interpr3.cxx

Go to the documentation of this file.

/* -*- 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 <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;
 
// 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 )
 
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::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
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
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
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;
}
 
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);
}
 
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
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
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;
    }
}
 
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(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 ) )
        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;
                    }
                }
                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());
            }
        }
    }
}
 
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 );
    }
}
 
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() );
}
 
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();
 
    KahanSum fSum    = 0.0;
    KahanSum 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(mrDoc, 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( mrContext, 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.get()/rValCount;
 
        if (nParamCount != 3)
        {
            sigma = (fSumSqr - fSum*fSum/rValCount).get()/(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;
    KahanSum fSum1    = 0.0;
    KahanSum fSumSqr1 = 0.0;
    KahanSum fSum2    = 0.0;
    KahanSum fSumSqr2 = 0.0;
    double fVal;
    SCSIZE i,j;
    for (i = 0; i < nC1; i++)
        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).get() / (fCount1-1.0) / fCount1;
        double fS2 = (fSumSqr2-fSum2*fSum2/fCount2).get() / (fCount2-1.0) / fCount2;
        if (fS1 + fS2 == 0.0)
        {
            PushNoValue();
            return false;
        }
        fT = std::abs(( fSum1/fCount1 - fSum2/fCount2 ).get())/sqrt(fS1+fS2);
        double c = fS1/(fS1+fS2);
    //  GetTDist is calculated via GetBetaDist and also works with non-integral
    // degrees of freedom. The result matches Excel
        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).get() / (fCount1 - 1.0);    // Variance
        double fS2 = (fSumSqr2 - fSum2*fSum2/fCount2).get() / (fCount2 - 1.0);
        fT = std::abs( fSum1.get()/fCount1 - fSum2.get()/fCount2 ) /
             sqrt( (fCount1-1.0)*fS1 + (fCount2-1.0)*fS2 ) *
             sqrt( fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2) );
        fF = fCount1 + fCount2 - 2;
    }
    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;
        KahanSum fSum1    = 0.0;
        KahanSum fSum2    = 0.0;
        KahanSum fSumSqrD = 0.0;
        double fVal1, fVal2;
        for (i = 0; i < nC1; i++)
            for (j = 0; j < nR1; j++)
            {
                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;
        }
        KahanSum fSumD = fSum1 - fSum2;
        double fDivider = ( fSumSqrD*fCount - fSumD*fSumD ).get();
        if ( fDivider == 0.0 )
        {
            PushError(FormulaError::DivisionByZero);
            return;
        }
        fT = std::abs(fSumD.get()) * 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;
    }
 
    auto aVal1 = pMat1->CollectKahan(sc::op::kOpSumAndSumSquare);
    auto aVal2 = pMat2->CollectKahan(sc::op::kOpSumAndSumSquare);
    double fCount1  = aVal1.mnCount;
    double fCount2  = aVal2.mnCount;
    KahanSum fSum1    = aVal1.maAccumulator[0];
    KahanSum fSumSqr1 = aVal1.maAccumulator[1];
    KahanSum fSum2    = aVal2.maAccumulator[0];
    KahanSum fSumSqr2 = aVal2.maAccumulator[1];
 
    if (fCount1 < 2.0 || fCount2 < 2.0)
    {
        PushNoValue();
        return;
    }
    double fS1 = (fSumSqr1-fSum1*fSum1/fCount1).get() / (fCount1-1.0);
    double fS2 = (fSumSqr2-fSum2*fSum2/fCount2).get() / (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;
    }
    KahanSum 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;
                    if (std::isinf(fTemp2))
                    {
                        PushError(FormulaError::NoConvergence);
                        return;
                    }
                    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.get(), fDF));
}
 
void ScInterpreter::ScKurt()
{
    KahanSum fSum;
    double fCount;
    std::vector<double> values;
    if ( !CalculateSkew(fSum, fCount, values) )
        return;
 
    // ODF 1.2 constraints: # of numbers >= 4
    if (fCount < 4.0)
    {
        // for interoperability with Excel
        PushError( FormulaError::DivisionByZero);
        return;
    }
 
    KahanSum vSum;
    double fMean = fSum.get() / fCount;
    for (double v : values)
        vSum += (v - fMean) * (v - fMean);
 
    double fStdDev = sqrt(vSum.get() / (fCount - 1.0));
    if (fStdDev == 0.0)
    {
        PushError( FormulaError::DivisionByZero);
        return;
    }
 
    KahanSum xpower4 = 0.0;
    for (double v : values)
    {
        double dx = (v - fMean) / fStdDev;
        xpower4 += (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.get() * k_l - k_t);
}
 
void ScInterpreter::ScHarMean()
{
    short nParamCount = GetByte();
    KahanSum 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(mrDoc, 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( mrContext, 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.get() );
    else
        PushError( nGlobalError);
}
 
void ScInterpreter::ScGeoMean()
{
    short nParamCount = GetByte();
    KahanSum 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(mrDoc, 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(mrContext, 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.get() / 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(KahanSum& fSum, double& fCount, std::vector<double>& values)
{
    short nParamCount = GetByte();
    if ( !MustHaveParamCountMin( nParamCount, 1 )  )
        return false;
 
    fSum   = 0.0;
    fCount = 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(mrDoc, 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( mrContext, 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 )
{
    KahanSum fSum;
    double fCount;
    std::vector<double> values;
    if (!CalculateSkew( fSum, fCount, values))
        return;
     // SKEW/SKEWP's constraints: they require at least three numbers
    if (fCount < 3.0)
    {
        // for interoperability with Excel
        PushError(FormulaError::DivisionByZero);
        return;
    }
 
    KahanSum vSum;
    double fMean = fSum.get() / fCount;
    for (double v : values)
        vSum += (v - fMean) * (v - fMean);
 
    double fStdDev = sqrt( vSum.get() / (bSkewp ? fCount : (fCount - 1.0)));
    if (fStdDev == 0)
    {
        PushIllegalArgument();
        return;
    }
 
    KahanSum xcube = 0.0;
    for (double v : values)
    {
        double dx = (v - fMean) / fStdDev;
        xcube += dx * dx * dx;
    }
 
    if (bSkewp)
        PushDouble( xcube.get() / fCount );
    else
        PushDouble( ((xcube.get() * 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)
        {
            // Note: neg fDiff seen with forum-mso-en4-719754.xlsx with
            // fPercentile of near 1 where approxFloor gave nIndex of nSize-1
            // resulting in a non-zero tiny negative fDiff.
            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 )
    {
        PushNoValue();
        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 )
    {
        PushNoValue();
        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;
 
    SCSIZE nCol = 0, nRow = 0;
    const auto aArray = GetRankNumberArray(nCol, nRow);
    const size_t nRankArraySize = aArray.size();
    if (nRankArraySize == 0 || nGlobalError != FormulaError::NONE)
    {
        PushNoValue();
        return;
    }
    assert(nRankArraySize == nCol * nRow);
 
    std::vector<SCSIZE> aRankArray;
    aRankArray.reserve(nRankArraySize);
    std::transform(aArray.begin(), aArray.end(), std::back_inserter(aRankArray),
            [bSmall](double f) {
                f = (bSmall ? rtl::math::approxFloor(f) : rtl::math::approxCeil(f));
                // Valid ranks are >= 1.
                if (f < 1.0 || !o3tl::convertsToAtMost(f, std::numeric_limits<SCSIZE>::max()))
                    return static_cast<SCSIZE>(0);
                return static_cast<SCSIZE>(f);
            });
 
    vector<double> aSortArray;
    GetNumberSequenceArray(1, aSortArray, false );
    const SCSIZE nSize = aSortArray.size();
    if (nSize == 0 || nGlobalError != FormulaError::NONE)
        PushNoValue();
    else if (nRankArraySize == 1)
    {
        const SCSIZE k = aRankArray[0];
        if (k < 1 || nSize < k)
        {
            if (!std::isfinite(aArray[0]))
                PushDouble(aArray[0]);  // propagates error
            else
                PushNoValue();
        }
        else
        {
            vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k);
            ::std::nth_element( aSortArray.begin(), iPos, aSortArray.end());
            PushDouble( *iPos);
        }
    }
    else
    {
        std::set<SCSIZE> aIndices;
        for (SCSIZE n : aRankArray)
        {
            if (1 <= n && n <= nSize)
                aIndices.insert(bSmall ? n-1 : nSize-n);
        }
        // We can spare sorting when the total number of ranks is small enough.
        // Find only the elements at given indices if, arbitrarily, the index size is
        // smaller than 1/3 of the haystack array's size; just sort it squarely, otherwise.
        if (aIndices.size() < nSize/3)
        {
            auto itBegin = aSortArray.begin();
            for (SCSIZE i : aIndices)
            {
                auto it = aSortArray.begin() + i;
                std::nth_element(itBegin, it, aSortArray.end());
                itBegin = ++it;
            }
        }
        else
            std::sort(aSortArray.begin(), aSortArray.end());
 
        std::vector<double> aResultArray;
        aResultArray.reserve(nRankArraySize);
        for (size_t i = 0; i < nRankArraySize; ++i)
        {
            const SCSIZE n = aRankArray[i];
            if (1 <= n && n <= nSize)
                aResultArray.push_back( aSortArray[bSmall ? n-1 : nSize-n]);
            else if (!std::isfinite( aArray[i]))
                aResultArray.push_back( aArray[i]);  // propagate error
            else
                aResultArray.push_back( CreateDoubleError( FormulaError::IllegalArgument));
        }
        ScMatrixRef pResult = GetNewMat(nCol, nRow, aResultArray);
        PushMatrix(pResult);
    }
}
 
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");
        KahanSum fSum = 0.0;
        for (SCSIZE i = nIndex; i < nSize-nIndex; i++)
            fSum += aSortArray[i];
        PushDouble(fSum.get()/static_cast<double>(nSize-2*nIndex));
    }
}
 
std::vector<double> ScInterpreter::GetRankNumberArray( SCSIZE& rCol, SCSIZE& rRow )
{
    std::vector<double> aArray;
    switch (GetStackType())
    {
        case svDouble:
            aArray.push_back(PopDouble());
            rCol = rRow = 1;
        break;
        case svSingleRef:
        {
            ScAddress aAdr;
            PopSingleRef(aAdr);
            ScRefCellValue aCell(mrDoc, aAdr);
            if (aCell.hasNumeric())
            {
                aArray.push_back(GetCellValue(aAdr, aCell));
                rCol = rRow = 1;
            }
        }
        break;
        case svDoubleRef:
        {
            ScRange aRange;
            PopDoubleRef(aRange, true);
            if (nGlobalError != FormulaError::NONE)
                break;
 
            // give up unless the start and end are in the same sheet
            if (aRange.aStart.Tab() != aRange.aEnd.Tab())
            {
                SetError(FormulaError::IllegalParameter);
                break;
            }
 
            // the range already is in order
            assert(aRange.aStart.Col() <= aRange.aEnd.Col());
            assert(aRange.aStart.Row() <= aRange.aEnd.Row());
            rCol = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
            rRow = aRange.aEnd.Row() - aRange.aStart.Row() + 1;
            aArray.reserve(rCol * rRow);
 
            FormulaError nErr = FormulaError::NONE;
            double fCellVal;
            ScValueIterator aValIter(mrContext, aRange, mnSubTotalFlags);
            if (aValIter.GetFirst(fCellVal, nErr))
            {
                do
                    aArray.push_back(fCellVal);
                while (aValIter.GetNext(fCellVal, nErr) && nErr == FormulaError::NONE);
            }
            // Note that SMALL() and LARGE() rank parameters (2nd) have
            // ParamClass::Value, so in array mode this is never hit and
            // argument was converted to matrix instead, but for normal
            // evaluation any non-numeric value including empty cell will
            // result in error anyway, so just clear and propagate an existing
            // error here already.
            if (aArray.size() != rCol * rRow)
            {
                aArray.clear();
                SetError(nErr);
            }
        }
        break;
        case svMatrix:
        case svExternalSingleRef:
        case svExternalDoubleRef:
        {
            ScMatrixRef pMat = GetMatrix();
            if (!pMat)
                break;
 
            const SCSIZE nCount = pMat->GetElementCount();
            aArray.reserve(nCount);
            // Do not propagate errors from matrix elements as global error.
            pMat->SetErrorInterpreter(nullptr);
            if (pMat->IsNumeric())
            {
                for (SCSIZE i = 0; i < nCount; ++i)
                    aArray.push_back(pMat->GetDouble(i));
            }
            else
            {
                for (SCSIZE i = 0; i < nCount; ++i)
                {
                    if (pMat->IsValue(i))
                        aArray.push_back( pMat->GetDouble(i));
                    else
                        aArray.push_back( CreateDoubleError( FormulaError::NoValue));
                }
            }
            pMat->GetDimensions(rCol, rRow);
        }
        break;
        default:
            PopError();
            SetError(FormulaError::IllegalParameter);
        break;
    }
    return aArray;
}
 
void ScInterpreter::GetNumberSequenceArray( sal_uInt8 nParamCount, vector<double>& rArray, bool bConvertTextInArray )
{
    ScAddress aAdr;
    ScRange aRange;
    const bool bIgnoreErrVal = bool(mnSubTotalFlags & SubtotalFlags::IgnoreErrVal);
    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(mrDoc, aAdr);
                if (bIgnoreErrVal && aCell.hasError())
                    ;   // nothing
                else 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( mrContext, aRange, mnSubTotalFlags );
                if (aValIter.GetFirst( fCellVal, nErr))
                {
                    if (bIgnoreErrVal)
                    {
                        if (nErr == FormulaError::NONE)
                            rArray.push_back( fCellVal);
                        while (aValIter.GetNext( fCellVal, nErr))
                        {
                            if (nErr == FormulaError::NONE)
                                rArray.push_back( fCellVal);
                        }
                    }
                    else
                    {
                        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())
                {
                    if (bIgnoreErrVal)
                    {
                        for (SCSIZE i = 0; i < nCount; ++i)
                        {
                            const double fVal = pMat->GetDouble(i);
                            if (nGlobalError == FormulaError::NONE)
                                rArray.push_back( fVal);
                            else
                                nGlobalError = FormulaError::NONE;
                        }
                    }
                    else
                    {
                        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 ) )
                        {
                            if (bIgnoreErrVal)
                            {
                                const double fVal = pMat->GetDouble(i);
                                if (nGlobalError == FormulaError::NONE)
                                    rArray.push_back( fVal);
                                else
                                    nGlobalError = FormulaError::NONE;
                            }
                            else
                                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
                                {
                                    if (!bIgnoreErrVal)
                                        rArray.push_back( CreateDoubleError( FormulaError::NoValue));
                                    // Propagate previous error if any, else
                                    // the current #VALUE! error, unless
                                    // ignoring error values.
                                    if (nErr != FormulaError::NONE)
                                        nGlobalError = nErr;
                                    else if (!bIgnoreErrVal)
                                        nGlobalError = FormulaError::NoValue;
                                    else
                                        nGlobalError = FormulaError::NONE;
                                }
                            }
                        }
                    }
                }
                else
                {
                    if (bIgnoreErrVal)
                    {
                        for (SCSIZE i = 0; i < nCount; ++i)
                        {
                            if (pMat->IsValue(i))
                            {
                                const double fVal = pMat->GetDouble(i);
                                if (nGlobalError == FormulaError::NONE)
                                    rArray.push_back( fVal);
                                else
                                    nGlobalError = FormulaError::NONE;
                            }
                        }
                    }
                    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<tools::Long>* pIndexOrder, bool bConvertTextInArray, bool bAllowEmptyArray )
{
    GetNumberSequenceArray( nParamCount, rSortArray, bConvertTextInArray );
    if (rSortArray.size() > MAX_COUNT_DOUBLE_FOR_SORT(mrDoc.GetSheetLimits()))
        SetError( FormulaError::MatrixSize);
    else if ( rSortArray.empty() )
    {
        if ( bAllowEmptyArray )
            return;
        SetError( FormulaError::NoValue);
    }
 
    if (nGlobalError == FormulaError::NONE)
        QuickSort( rSortArray, pIndexOrder);
}
 
static void lcl_QuickSort( tools::Long nLo, tools::Long nHi, vector<double>& rSortArray, vector<tools::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;
    }
 
    tools::Long ni = nLo;
    tools::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<tools::Long>* pIndexOrder )
{
    tools::Long n = static_cast<tools::Long>(rSortArray.size());
 
    if (pIndexOrder)
    {
        pIndexOrder->clear();
        pIndexOrder->reserve(n);
        for (tools::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 ] )
            PushError( FormulaError::NotAvailable);
        else
        {
            double fLastPos = 0;
            double fFirstPos = -1.0;
            bool bFinished = false;
            SCSIZE i;
            for (i = 0; i < nSize && !bFinished; 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 (fFirstPos <= 0)
            {
                PushError( FormulaError::NotAvailable);
            }
            else 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;
    KahanSum 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(mrDoc, 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( mrContext, 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.get() / rValCount;
    sp = SaveSP;
    rVal = 0.0;
    nParam = nParamCount;
    nRefInList = 0;
    while (nParam-- > 0)
    {
        switch (GetStackType())
        {
            case svDouble :
                rVal += std::abs(GetDouble() - nMiddle);
                break;
            case svSingleRef :
            {
                PopSingleRef( aAdr );
                ScRefCellValue aCell(mrDoc, aAdr);
                if (aCell.hasNumeric())
                    rVal += std::abs(GetCellValue(aAdr, aCell) - nMiddle);
            }
            break;
            case svDoubleRef :
            case svRefList :
            {
                FormulaError nErr = FormulaError::NONE;
                double nCellVal;
                PopDoubleRef( aRange, nParam, nRefInList);
                ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
                if (aValIter.GetFirst(nCellVal, nErr))
                {
                    rVal += std::abs(nCellVal - nMiddle);
                    while (aValIter.GetNext(nCellVal, nErr))
                         rVal += std::abs(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 += std::abs(pMat->GetDouble(nElem) - nMiddle);
                        }
                    }
                    else
                    {
                        for (SCSIZE nElem = 0; nElem < nCount; nElem++)
                        {
                            if (!pMat->IsStringOrEmpty(nElem))
                                rVal += std::abs(pMat->GetDouble(nElem) - nMiddle);
                        }
                    }
                }
            }
            break;
            default : SetError(FormulaError::IllegalParameter); break;
        }
    }
    PushDouble(rVal.get() / 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)
        std::swap( fLo, fUp );
    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
        {
            KahanSum fSum = 0.0;
            KahanSum 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 || std::abs((fSum -1.0).get()) > 1.0E-7)
                PushNoValue();
            else
                PushDouble(fRes.get());
        }
    }
}
 
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;
    KahanSum fSumX          = 0.0;
    KahanSum 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))
            {
                fSumX += pMat1->GetDouble(i,j);
                fSumY += pMat2->GetDouble(i,j);
                fCount++;
            }
        }
    }
    if (fCount < (_bStexy ? 3.0 : (_bSample ? 2.0 : 1.0)))
        PushNoValue();
    else
    {
        KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
        KahanSum fSumSqrDeltaX    = 0.0; // sum of (ValX-MeanX)^2
        KahanSum fSumSqrDeltaY    = 0.0; // sum of (ValY-MeanY)^2
        const double fMeanX = fSumX.get() / fCount;
        const double fMeanY = fSumY.get() / 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 )
        {
            // tdf#94962 - Values below the numerical limit lead to unexpected results
            if (fSumSqrDeltaX < ::std::numeric_limits<double>::min()
                || (!_bStexy && fSumSqrDeltaY < ::std::numeric_limits<double>::min()))
                PushError( FormulaError::DivisionByZero);
            else if ( _bStexy )
                PushDouble( sqrt( ( fSumSqrDeltaY - fSumDeltaXDeltaY *
                            fSumDeltaXDeltaY / fSumSqrDeltaX ).get() / (fCount-2)));
            else
                PushDouble( fSumDeltaXDeltaY.get() / sqrt( fSumSqrDeltaX.get() * fSumSqrDeltaY.get() ));
        }
        else
        {
            if ( _bSample )
                PushDouble( fSumDeltaXDeltaY.get() / ( fCount - 1 ) );
            else
                PushDouble( fSumDeltaXDeltaY.get() / fCount);
        }
    }
}
 
void ScInterpreter::ScRSQ()
{
    // Same as ScPearson()*ScPearson()
    ScPearson();
    if (nGlobalError != FormulaError::NONE)
        return;
 
    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;
    KahanSum fSumX          = 0.0;
    KahanSum 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))
            {
                fSumX += pMat1->GetDouble(i,j);
                fSumY += pMat2->GetDouble(i,j);
                fCount++;
            }
        }
    }
    if (fCount < 1.0)
        PushNoValue();
    else
    {
        KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
        KahanSum fSumSqrDeltaX    = 0.0; // sum of (ValX-MeanX)^2
        double fMeanX = fSumX.get() / fCount;
        double fMeanY = fSumY.get() / 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.get() / fSumSqrDeltaX.get());
            else
                PushDouble( fMeanY - fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get() * 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;
    KahanSum fSumX          = 0.0;
    KahanSum 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))
            {
                fSumX += pMat1->GetDouble(i,j);
                fSumY += pMat2->GetDouble(i,j);
                fCount++;
            }
        }
    }
    if (fCount < 1.0)
        PushNoValue();
    else
    {
        KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
        KahanSum fSumSqrDeltaX    = 0.0; // sum of (ValX-MeanX)^2
        double fMeanX = fSumX.get() / fCount;
        double fMeanY = fSumY.get() / 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.get() / fSumSqrDeltaX.get() * (fVal - fMeanX));
    }
}
 
static void lcl_roundUpNearestPow2(SCSIZE& nNum, SCSIZE& nNumBits)
{
    // Find the least power of 2 that is less than or equal to nNum.
    SCSIZE nPow2(1);
    nNumBits = std::numeric_limits<SCSIZE>::digits;
    nPow2 <<= (nNumBits - 1);
    while (nPow2 >= 1)
    {
        if (nNum & nPow2)
            break;
 
        --nNumBits;
        nPow2 >>= 1;
    }
 
    if (nPow2 != nNum)
        nNum = nPow2 ? (nPow2 << 1) : 1;
    else
        --nNumBits;
}
 
static SCSIZE lcl_bitReverse(SCSIZE nIn, SCSIZE nBound)
{
    SCSIZE nOut = 0;
    for (SCSIZE nMask = 1; nMask < nBound; nMask <<= 1)
    {
        nOut <<= 1;
 
        if (nIn & nMask)
            nOut |= 1;
    }
 
    return nOut;
}
 
namespace {
 
// Computes and stores twiddle factors for computing DFT later.
struct ScTwiddleFactors
{
    ScTwiddleFactors(SCSIZE nN, bool bInverse) :
        mfWReal(nN),
        mfWImag(nN),
        mnN(nN),
        mbInverse(bInverse)
    {}
 
    void Compute();
 
    void Conjugate()
    {
        mbInverse = !mbInverse;
        for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
            mfWImag[nIdx] = -mfWImag[nIdx];
    }
 
    std::vector<double> mfWReal;
    std::vector<double> mfWImag;
    SCSIZE mnN;
    bool mbInverse;
};
 
}
 
void ScTwiddleFactors::Compute()
{
    mfWReal.resize(mnN);
    mfWImag.resize(mnN);
 
    double nW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnN);
 
    if (mnN == 1)
    {
        mfWReal[0] = 1.0;
        mfWImag[0] = 0.0;
    }
    else if (mnN == 2)
    {
        mfWReal[0] = 1;
        mfWImag[0] = 0;
 
        mfWReal[1] = -1;
        mfWImag[1] = 0;
    }
    else if (mnN == 4)
    {
        mfWReal[0] = 1;
        mfWImag[0] = 0;
 
        mfWReal[1] = 0;
        mfWImag[1] = (mbInverse ? 1.0 : -1.0);
 
        mfWReal[2] = -1;
        mfWImag[2] = 0;
 
        mfWReal[3] = 0;
        mfWImag[3] = (mbInverse ? -1.0 : 1.0);
    }
    else if ((mnN % 4) == 0)
    {
        const SCSIZE nQSize = mnN >> 2;
        // Compute cos of the start quadrant.
        // This is the first quadrant if mbInverse == true, else it is the fourth quadrant.
        for (SCSIZE nIdx = 0; nIdx <= nQSize; ++nIdx)
            mfWReal[nIdx] = cos(nW*static_cast<double>(nIdx));
 
        if (mbInverse)
        {
            const SCSIZE nQ1End = nQSize;
            // First quadrant
            for (SCSIZE nIdx = 0; nIdx <= nQ1End; ++nIdx)
                mfWImag[nIdx] = mfWReal[nQ1End-nIdx];
 
            // Second quadrant
            const SCSIZE nQ2End = nQ1End << 1;
            for (SCSIZE nIdx = nQ1End+1; nIdx <= nQ2End; ++nIdx)
            {
                mfWReal[nIdx] = -mfWReal[nQ2End - nIdx];
                mfWImag[nIdx] =  mfWImag[nQ2End - nIdx];
            }
 
            // Third quadrant
            const SCSIZE nQ3End = nQ2End + nQ1End;
            for (SCSIZE nIdx = nQ2End+1; nIdx <= nQ3End; ++nIdx)
            {
                mfWReal[nIdx] = -mfWReal[nIdx - nQ2End];
                mfWImag[nIdx] = -mfWImag[nIdx - nQ2End];
            }
 
            // Fourth Quadrant
            for (SCSIZE nIdx = nQ3End+1; nIdx < mnN; ++nIdx)
            {
                mfWReal[nIdx] =  mfWReal[mnN - nIdx];
                mfWImag[nIdx] = -mfWImag[mnN - nIdx];
            }
        }
        else
        {
            const SCSIZE nQ4End = nQSize;
            const SCSIZE nQ3End = nQSize << 1;
            const SCSIZE nQ2End = nQ3End + nQSize;
 
            // Fourth quadrant.
            for (SCSIZE nIdx = 0; nIdx <= nQ4End; ++nIdx)
                mfWImag[nIdx] = -mfWReal[nQ4End - nIdx];
 
            // Third quadrant.
            for (SCSIZE nIdx = nQ4End+1; nIdx <= nQ3End; ++nIdx)
            {
                mfWReal[nIdx] = -mfWReal[nQ3End - nIdx];
                mfWImag[nIdx] =  mfWImag[nQ3End - nIdx];
            }
 
            // Second quadrant.
            for (SCSIZE nIdx = nQ3End+1; nIdx <= nQ2End; ++nIdx)
            {
                mfWReal[nIdx] = -mfWReal[nIdx - nQ3End];
                mfWImag[nIdx] = -mfWImag[nIdx - nQ3End];
            }
 
            // First quadrant.
            for (SCSIZE nIdx = nQ2End+1; nIdx < mnN; ++nIdx)
            {
                mfWReal[nIdx] =  mfWReal[mnN - nIdx];
                mfWImag[nIdx] = -mfWImag[mnN - nIdx];
            }
        }
    }
    else
    {
        for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
        {
            double fAngle = nW*static_cast<double>(nIdx);
            mfWReal[nIdx] = cos(fAngle);
            mfWImag[nIdx] = sin(fAngle);
        }
    }
}
 
namespace {
 
// A radix-2 decimation in time FFT algorithm for complex valued input.
class ScComplexFFT2
{
public:
    // rfArray.size() would always be even and a power of 2. (asserted in prepare())
    // rfArray's first half contains the real parts and the later half contains the imaginary parts.
    ScComplexFFT2(std::vector<double>& raArray, bool bInverse, bool bPolar, double fMinMag,
                  ScTwiddleFactors& rTF, bool bSubSampleTFs = false, bool bDisableNormalize = false) :
        mrArray(raArray),
        mfWReal(rTF.mfWReal),
        mfWImag(rTF.mfWImag),
        mnPoints(raArray.size()/2),
        mnStages(0),
        mfMinMag(fMinMag),
        mbInverse(bInverse),
        mbPolar(bPolar),
        mbDisableNormalize(bDisableNormalize),
        mbSubSampleTFs(bSubSampleTFs)
    {}
 
    void Compute();
 
private:
 
    void prepare();
 
    double getReal(SCSIZE nIdx)
    {
        return mrArray[nIdx];
    }
 
    void setReal(double fVal, SCSIZE nIdx)
    {
        mrArray[nIdx] = fVal;
    }
 
    double getImag(SCSIZE nIdx)
    {
        return mrArray[mnPoints + nIdx];
    }
 
    void setImag(double fVal, SCSIZE nIdx)
    {
        mrArray[mnPoints + nIdx] = fVal;
    }
 
    SCSIZE getTFactorIndex(SCSIZE nPtIndex, SCSIZE nTfIdxScaleBits)
    {
        return ( ( nPtIndex << nTfIdxScaleBits ) & ( mnPoints - 1 ) ); // (x & (N-1)) is same as (x % N) but faster.
    }
 
    void computeFly(SCSIZE nTopIdx, SCSIZE nBottomIdx, SCSIZE nWIdx1, SCSIZE nWIdx2)
    {
        if (mbSubSampleTFs)
        {
            nWIdx1 <<= 1;
            nWIdx2 <<= 1;
        }
 
        const double x1r = getReal(nTopIdx);
        const double x2r = getReal(nBottomIdx);
 
        const double& w1r = mfWReal[nWIdx1];
        const double& w1i = mfWImag[nWIdx1];
 
        const double& w2r = mfWReal[nWIdx2];
        const double& w2i = mfWImag[nWIdx2];
 
        const double x1i = getImag(nTopIdx);
        const double x2i = getImag(nBottomIdx);
 
        setReal(x1r + x2r*w1r - x2i*w1i, nTopIdx);
        setImag(x1i + x2i*w1r + x2r*w1i, nTopIdx);
 
        setReal(x1r + x2r*w2r - x2i*w2i, nBottomIdx);
        setImag(x1i + x2i*w2r + x2r*w2i, nBottomIdx);
    }
 
    std::vector<double>& mrArray;
    std::vector<double>& mfWReal;
    std::vector<double>& mfWImag;
    SCSIZE mnPoints;
    SCSIZE mnStages;
    double mfMinMag;
    bool mbInverse:1;
    bool mbPolar:1;
    bool mbDisableNormalize:1;
    bool mbSubSampleTFs:1;
};
 
}
 
void ScComplexFFT2::prepare()
{
    SCSIZE nPoints = mnPoints;
    lcl_roundUpNearestPow2(nPoints, mnStages);
    assert(nPoints == mnPoints);
 
    // Reorder array by bit-reversed indices.
    for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
    {
        SCSIZE nRevIdx = lcl_bitReverse(nIdx, mnPoints);
        if (nIdx < nRevIdx)
        {
            double fTmp = getReal(nIdx);
            setReal(getReal(nRevIdx), nIdx);
            setReal(fTmp, nRevIdx);
 
            fTmp = getImag(nIdx);
            setImag(getImag(nRevIdx), nIdx);
            setImag(fTmp, nRevIdx);
        }
    }
}
 
static void lcl_normalize(std::vector<double>& rCmplxArray, bool bScaleOnlyReal)
{
    const SCSIZE nPoints = rCmplxArray.size()/2;
    const double fScale = 1.0/static_cast<double>(nPoints);
 
    // Scale the real part
    for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
        rCmplxArray[nIdx] *= fScale;
 
    if (!bScaleOnlyReal)
    {
        const SCSIZE nLen = nPoints*2;
        for (SCSIZE nIdx = nPoints; nIdx < nLen; ++nIdx)
            rCmplxArray[nIdx] *= fScale;
    }
}
 
static void lcl_convertToPolar(std::vector<double>& rCmplxArray, double fMinMag)
{
    const SCSIZE nPoints = rCmplxArray.size()/2;
    double fMag, fPhase, fR, fI;
    for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
    {
        fR = rCmplxArray[nIdx];
        fI = rCmplxArray[nPoints+nIdx];
        fMag = std::hypot(fR, fI);
        if (fMag < fMinMag)
        {
            fMag = 0.0;
            fPhase = 0.0;
        }
        else
        {
            fPhase = atan2(fI, fR);
        }
 
        rCmplxArray[nIdx] = fMag;
        rCmplxArray[nPoints+nIdx] = fPhase;
    }
}
 
void ScComplexFFT2::Compute()
{
    prepare();
 
    const SCSIZE nFliesInStage = mnPoints/2;
    for (SCSIZE nStage = 0; nStage < mnStages; ++nStage)
    {
        const SCSIZE nTFIdxScaleBits = mnStages - nStage - 1;  // Twiddle factor index's scale factor in bits.
        const SCSIZE nFliesInGroup = SCSIZE(1) << nStage;
        const SCSIZE nGroups = nFliesInStage/nFliesInGroup;
        const SCSIZE nFlyWidth = nFliesInGroup;
        for (SCSIZE nGroup = 0, nFlyTopIdx = 0; nGroup < nGroups; ++nGroup)
        {
            for (SCSIZE nFly = 0; nFly < nFliesInGroup; ++nFly, ++nFlyTopIdx)
            {
                SCSIZE nFlyBottomIdx = nFlyTopIdx + nFlyWidth;
                SCSIZE nWIdx1 = getTFactorIndex(nFlyTopIdx, nTFIdxScaleBits);
                SCSIZE nWIdx2 = getTFactorIndex(nFlyBottomIdx, nTFIdxScaleBits);
 
                computeFly(nFlyTopIdx, nFlyBottomIdx, nWIdx1, nWIdx2);
            }
 
            nFlyTopIdx += nFlyWidth;
        }
    }
 
    if (mbPolar)
        lcl_convertToPolar(mrArray, mfMinMag);
 
    // Normalize after converting to polar, so we have a chance to
    // save O(mnPoints) flops.
    if (mbInverse && !mbDisableNormalize)
        lcl_normalize(mrArray, mbPolar);
}
 
namespace {
 
// Bluestein's algorithm or chirp z-transform algorithm that can be used to
// compute DFT of a complex valued input of any length N in O(N lgN) time.
class ScComplexBluesteinFFT
{
public:
 
    ScComplexBluesteinFFT(std::vector<double>& rArray, bool bReal, bool bInverse,
                          bool bPolar, double fMinMag, bool bDisableNormalize = false) :
        mrArray(rArray),
        mnPoints(rArray.size()/2), // rArray should have space for imaginary parts even if real input.
        mfMinMag(fMinMag),
        mbReal(bReal),
        mbInverse(bInverse),
        mbPolar(bPolar),
        mbDisableNormalize(bDisableNormalize)
    {}
 
    void Compute();
 
private:
    std::vector<double>& mrArray;
    const SCSIZE mnPoints;
    double mfMinMag;
    bool mbReal:1;
    bool mbInverse:1;
    bool mbPolar:1;
    bool mbDisableNormalize:1;
};
 
}
 
void ScComplexBluesteinFFT::Compute()
{
    std::vector<double> aRealScalars(mnPoints);
    std::vector<double> aImagScalars(mnPoints);
    double fW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnPoints);
    for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
    {
        double fAngle = 0.5*fW*static_cast<double>(nIdx*nIdx);
        aRealScalars[nIdx] = cos(fAngle);
        aImagScalars[nIdx] = sin(fAngle);
    }
 
    SCSIZE nMinSize = mnPoints*2 - 1;
    SCSIZE nExtendedLength = nMinSize, nTmp = 0;
    lcl_roundUpNearestPow2(nExtendedLength, nTmp);
    std::vector<double> aASignal(nExtendedLength*2); // complex valued
    std::vector<double> aBSignal(nExtendedLength*2); // complex valued
 
    double fReal, fImag;
    for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
    {
        // Real part of A signal.
        aASignal[nIdx] = mrArray[nIdx]*aRealScalars[nIdx] + (mbReal ? 0.0 : -mrArray[mnPoints+nIdx]*aImagScalars[nIdx]);
        // Imaginary part of A signal.
        aASignal[nExtendedLength + nIdx] = mrArray[nIdx]*aImagScalars[nIdx] + (mbReal ? 0.0 : mrArray[mnPoints+nIdx]*aRealScalars[nIdx]);
 
        // Real part of B signal.
        aBSignal[nIdx] = fReal = aRealScalars[nIdx];
        // Imaginary part of B signal.
        aBSignal[nExtendedLength + nIdx] = fImag = -aImagScalars[nIdx]; // negative sign because B signal is the conjugation of the scalars.
 
        if (nIdx)
        {
            // B signal needs a mirror of its part in 0 < n < mnPoints at the tail end.
            aBSignal[nExtendedLength - nIdx] = fReal;
            aBSignal[(nExtendedLength<<1) - nIdx] = fImag;
        }
    }
 
    {
        ScTwiddleFactors aTF(nExtendedLength, false /*not inverse*/);
        aTF.Compute();
 
        // Do complex-FFT2 of both A and B signal.
        ScComplexFFT2 aFFT2A(aASignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
                             aTF, false /*no subsample*/, true /* disable normalize */);
        aFFT2A.Compute();
 
        ScComplexFFT2 aFFT2B(aBSignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
                             aTF, false /*no subsample*/, true /* disable normalize */);
        aFFT2B.Compute();
 
        double fAR, fAI, fBR, fBI;
        for (SCSIZE nIdx = 0; nIdx < nExtendedLength; ++nIdx)
        {
            fAR = aASignal[nIdx];
            fAI = aASignal[nExtendedLength + nIdx];
            fBR = aBSignal[nIdx];
            fBI = aBSignal[nExtendedLength + nIdx];
 
            // Do point-wise product inplace in A signal.
            aASignal[nIdx] = fAR*fBR - fAI*fBI;
            aASignal[nExtendedLength + nIdx] = fAR*fBI + fAI*fBR;
        }
 
        // Do complex-inverse-FFT2 of aASignal.
        aTF.Conjugate();
        ScComplexFFT2 aFFT2AI(aASignal, true /*inverse*/, false /*no polar*/, 0.0 /* no clipping */, aTF); // Need normalization here.
        aFFT2AI.Compute();
    }
 
    // Point-wise multiply with scalars.
    for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
    {
        fReal = aASignal[nIdx];
        fImag = aASignal[nExtendedLength + nIdx];
        mrArray[nIdx] = fReal*aRealScalars[nIdx] - fImag*aImagScalars[nIdx]; // no conjugation needed here.
        mrArray[mnPoints + nIdx] = fReal*aImagScalars[nIdx] + fImag*aRealScalars[nIdx];
    }
 
    // Normalize/Polar operations
    if (mbPolar)
        lcl_convertToPolar(mrArray, mfMinMag);
 
    // Normalize after converting to polar, so we have a chance to
    // save O(mnPoints) flops.
    if (mbInverse && !mbDisableNormalize)
        lcl_normalize(mrArray, mbPolar);
}
 
namespace {
 
// Computes DFT of an even length(N) real-valued input by using a
// ScComplexFFT2 if N == 2^k for some k or else by using a ScComplexBluesteinFFT
// with a complex valued input of length = N/2.
class ScRealFFT
{
public:
 
    ScRealFFT(std::vector<double>& rInArray, std::vector<double>& rOutArray, bool bInverse,
              bool bPolar, double fMinMag) :
        mrInArray(rInArray),
        mrOutArray(rOutArray),
        mfMinMag(fMinMag),
        mbInverse(bInverse),
        mbPolar(bPolar)
    {}
 
    void Compute();
 
private:
    std::vector<double>& mrInArray;
    std::vector<double>& mrOutArray;
    double mfMinMag;
    bool mbInverse:1;
    bool mbPolar:1;
};
 
}
 
void ScRealFFT::Compute()
{
    // input length has to be even to do this optimization.
    assert(mrInArray.size() % 2 == 0);
    assert(mrInArray.size()*2 == mrOutArray.size());
    // nN is the number of points in the complex-fft input
    // which will be half of the number of points in real array.
    const SCSIZE nN = mrInArray.size()/2;
    if (nN == 0)
    {
        mrOutArray[0] = mrInArray[0];
        mrOutArray[1] = 0.0;
        return;
    }
 
    // work array should be the same length as mrInArray
    std::vector<double> aWorkArray(nN*2);
    for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
    {
        SCSIZE nDoubleIdx = 2*nIdx;
        // Use even elements as real part
        aWorkArray[nIdx] = mrInArray[nDoubleIdx];
        // and odd elements as imaginary part of the contrived complex sequence.
        aWorkArray[nN+nIdx] = mrInArray[nDoubleIdx+1];
    }
 
    ScTwiddleFactors aTFs(nN*2, mbInverse);
    aTFs.Compute();
    SCSIZE nNextPow2 = nN, nTmp = 0;
    lcl_roundUpNearestPow2(nNextPow2, nTmp);
 
    if (nNextPow2 == nN)
    {
        ScComplexFFT2 aFFT2(aWorkArray, mbInverse, false /*disable polar*/, 0.0 /* no clipping */,
                            aTFs, true /*subsample tf*/, true /*disable normalize*/);
        aFFT2.Compute();
    }
    else
    {
        ScComplexBluesteinFFT aFFT(aWorkArray, false /*complex input*/, mbInverse, false /*disable polar*/,
                                   0.0 /* no clipping */, true /*disable normalize*/);
        aFFT.Compute();
    }
 
    // Post process aWorkArray to populate mrOutArray
 
    const SCSIZE nTwoN = 2*nN, nThreeN = 3*nN;
    double fY1R, fY2R, fY1I, fY2I, fResR, fResI, fWR, fWI;
    for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
    {
        const SCSIZE nIdxRev = nIdx ? (nN - nIdx) : 0;
        fY1R = aWorkArray[nIdx];
        fY2R = aWorkArray[nIdxRev];
        fY1I = aWorkArray[nN + nIdx];
        fY2I = aWorkArray[nN + nIdxRev];
        fWR  = aTFs.mfWReal[nIdx];
        fWI  = aTFs.mfWImag[nIdx];
 
        // mrOutArray has length = 4*nN
        // Real part of the final output (only half of the symmetry around Nyquist frequency)
        // Fills the first quarter.
        mrOutArray[nIdx] = fResR = 0.5*(
            fY1R + fY2R +
            fWR * (fY1I + fY2I) +
            fWI * (fY1R - fY2R) );
        // Imaginary part of the final output (only half of the symmetry around Nyquist frequency)
        // Fills the third quarter.
        mrOutArray[nTwoN + nIdx] = fResI = 0.5*(
            fY1I - fY2I +
            fWI * (fY1I + fY2I) -
            fWR * (fY1R - fY2R) );
 
        // Fill the missing 2 quarters using symmetry argument.
        if (nIdx)
        {
            // Fills the 2nd quarter.
            mrOutArray[nN + nIdxRev] = fResR;
            // Fills the 4th quarter.
            mrOutArray[nThreeN + nIdxRev] = -fResI;
        }
        else
        {
            mrOutArray[nN] = fY1R - fY1I;
            mrOutArray[nThreeN] = 0.0;
        }
    }
 
    // Normalize/Polar operations
    if (mbPolar)
        lcl_convertToPolar(mrOutArray, mfMinMag);
 
    // Normalize after converting to polar, so we have a chance to
    // save O(mnPoints) flops.
    if (mbInverse)
        lcl_normalize(mrOutArray, mbPolar);
}
 
using ScMatrixGenerator = ScMatrixRef(SCSIZE, SCSIZE, std::vector<double>&);
 
namespace {
 
// Generic FFT class that decides which FFT implementation to use.
class ScFFT
{
public:
 
    ScFFT(ScMatrixRef& pMat, bool bReal, bool bInverse, bool bPolar, double fMinMag) :
        mpInputMat(pMat),
        mfMinMag(fMinMag),
        mbReal(bReal),
        mbInverse(bInverse),
        mbPolar(bPolar)
    {}
 
    ScMatrixRef Compute(const std::function<ScMatrixGenerator>& rMatGenFunc);
 
private:
    ScMatrixRef& mpInputMat;
    double mfMinMag;
    bool mbReal:1;
    bool mbInverse:1;
    bool mbPolar:1;
};
 
}
 
ScMatrixRef ScFFT::Compute(const std::function<ScMatrixGenerator>& rMatGenFunc)
{
    std::vector<double> aArray;
    mpInputMat->GetDoubleArray(aArray);
    SCSIZE nPoints = mbReal ? aArray.size() : (aArray.size()/2);
    if (nPoints == 1)
    {
        std::vector<double> aOutArray(2);
        aOutArray[0] = aArray[0];
        aOutArray[1] = mbReal ? 0.0 : aArray[1];
        if (mbPolar)
            lcl_convertToPolar(aOutArray, mfMinMag);
        return rMatGenFunc(2, 1, aOutArray);
    }
 
    if (mbReal && (nPoints % 2) == 0)
    {
        std::vector<double> aOutArray(nPoints*2);
        ScRealFFT aFFT(aArray, aOutArray, mbInverse, mbPolar, mfMinMag);
        aFFT.Compute();
        return rMatGenFunc(2, nPoints, aOutArray);
    }
 
    SCSIZE nNextPow2 = nPoints, nTmp = 0;
    lcl_roundUpNearestPow2(nNextPow2, nTmp);
    if (nNextPow2 == nPoints && !mbReal)
    {
        ScTwiddleFactors aTF(nPoints, mbInverse);
        aTF.Compute();
        ScComplexFFT2 aFFT2(aArray, mbInverse, mbPolar, mfMinMag, aTF);
        aFFT2.Compute();
        return rMatGenFunc(2, nPoints, aArray);
    }
 
    if (mbReal)
        aArray.resize(nPoints*2, 0.0);
    ScComplexBluesteinFFT aFFT(aArray, mbReal, mbInverse, mbPolar, mfMinMag);
    aFFT.Compute();
    return rMatGenFunc(2, nPoints, aArray);
}
 
void ScInterpreter::ScFourier()
{
    sal_uInt8 nParamCount = GetByte();
    if ( !MustHaveParamCount( nParamCount, 2, 5 ) )
        return;
 
    bool bInverse = false;
    bool bPolar = false;
    double fMinMag = 0.0;
 
    if (nParamCount == 5)
    {
        if (IsMissing())
            Pop();
        else
            fMinMag = GetDouble();
    }
 
    if (nParamCount >= 4)
    {
        if (IsMissing())
            Pop();
        else
            bPolar = GetBool();
    }
 
    if (nParamCount >= 3)
    {
        if (IsMissing())
            Pop();
        else
            bInverse = GetBool();
    }
 
    bool bGroupedByColumn = GetBool();
 
    ScMatrixRef pInputMat = GetMatrix();
    if (!pInputMat)
    {
        PushIllegalParameter();
        return;
    }
 
    SCSIZE nC, nR;
    pInputMat->GetDimensions(nC, nR);
 
    if ((bGroupedByColumn && nC > 2) || (!bGroupedByColumn && nR > 2))
    {
        // There can be no more than 2 columns (real, imaginary) if data grouped by columns.
        // and no more than 2 rows if data is grouped by rows.
        PushIllegalArgument();
        return;
    }
 
    if (!pInputMat->IsNumeric())
    {
        PushNoValue();
        return;
    }
 
    bool bRealInput = true;
    if (!bGroupedByColumn)
    {
        pInputMat->MatTrans(*pInputMat);
        bRealInput = (nR == 1);
    }
    else
    {
        bRealInput = (nC == 1);
    }
 
    ScFFT aFFT(pInputMat, bRealInput, bInverse, bPolar, fMinMag);
    std::function<ScMatrixGenerator> aFunc = [this](SCSIZE nCol, SCSIZE nRow, std::vector<double>& rVec) -> ScMatrixRef
    {
        return this->GetNewMat(nCol, nRow, rVec);
    };
    ScMatrixRef pOut = aFFT.Compute(aFunc);
    PushMatrix(pOut);
}
 
/* vim:set shiftwidth=4 softtabstop=4 expandtab: */