diff options
Diffstat (limited to 'sc/source/core/tool/interpr3.cxx')
| -rw-r--r-- | sc/source/core/tool/interpr3.cxx | 5579 |
1 files changed, 5579 insertions, 0 deletions
diff --git a/sc/source/core/tool/interpr3.cxx b/sc/source/core/tool/interpr3.cxx new file mode 100644 index 000000000..1fa4500a9 --- /dev/null +++ b/sc/source/core/tool/interpr3.cxx @@ -0,0 +1,5579 @@ +/* -*- 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; + +/// Two columns of data should be sortable with GetSortArray() and QuickSort() +// This is an arbitrary limit. +static size_t MAX_COUNT_DOUBLE_FOR_SORT(const ScSheetLimits& rSheetLimits) +{ + return rSheetLimits.GetMaxRowCount() * 2; +} + +const double ScInterpreter::fMaxGammaArgument = 171.624376956302; // found experimental +const double fMachEps = ::std::numeric_limits<double>::epsilon(); + +namespace { + +class ScDistFunc +{ +public: + virtual double GetValue(double x) const = 0; + +protected: + ~ScDistFunc() {} +}; + +} + +// iteration for inverse distributions + +//template< class T > double lcl_IterateInverse( const T& rFunction, double x0, double x1, bool& rConvError ) + +/** u*w<0.0 fails for values near zero */ +static bool lcl_HasChangeOfSign( double u, double w ) +{ + return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0); +} + +static double lcl_IterateInverse( const ScDistFunc& rFunction, double fAx, double fBx, bool& rConvError ) +{ + rConvError = false; + const double fYEps = 1.0E-307; + const double fXEps = ::std::numeric_limits<double>::epsilon(); + + OSL_ENSURE(fAx<fBx, "IterateInverse: wrong interval"); + + // find enclosing interval + + KahanSum fkAx = fAx; + KahanSum fkBx = fBx; + double fAy = rFunction.GetValue(fAx); + double fBy = rFunction.GetValue(fBx); + KahanSum fTemp; + unsigned short nCount; + for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy); nCount++) + { + if (std::abs(fAy) <= std::abs(fBy)) + { + fTemp = fkAx; + fkAx += (fkAx - fkBx) * 2.0; + if (fkAx < 0.0) + fkAx = 0.0; + fkBx = fTemp; + fBy = fAy; + fAy = rFunction.GetValue(fkAx.get()); + } + else + { + fTemp = fkBx; + fkBx += (fkBx - fkAx) * 2.0; + fkAx = fTemp; + fAy = fBy; + fBy = rFunction.GetValue(fkBx.get()); + } + } + + fAx = fkAx.get(); + fBx = fkBx.get(); + if (fAy == 0.0) + return fAx; + if (fBy == 0.0) + return fBx; + if (!lcl_HasChangeOfSign( fAy, fBy)) + { + rConvError = true; + return 0.0; + } + // inverse quadric interpolation with additional brackets + // set three points + double fPx = fAx; + double fPy = fAy; + double fQx = fBx; + double fQy = fBy; + double fRx = fAx; + double fRy = fAy; + double fSx = 0.5 * (fAx + fBx); // potential next point + bool bHasToInterpolate = true; + nCount = 0; + while ( nCount < 500 && std::abs(fRy) > fYEps && + (fBx-fAx) > ::std::max( std::abs(fAx), std::abs(fBx)) * fXEps ) + { + if (bHasToInterpolate) + { + if (fPy!=fQy && fQy!=fRy && fRy!=fPy) + { + fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy) + + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy) + + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy); + bHasToInterpolate = (fAx < fSx) && (fSx < fBx); // inside the brackets? + } + else + bHasToInterpolate = false; + } + if(!bHasToInterpolate) + { + fSx = 0.5 * (fAx + fBx); + // reset points + fQx = fBx; fQy = fBy; + bHasToInterpolate = true; + } + // shift points for next interpolation + fPx = fQx; fQx = fRx; fRx = fSx; + fPy = fQy; fQy = fRy; fRy = rFunction.GetValue(fSx); + // update brackets + if (lcl_HasChangeOfSign( fAy, fRy)) + { + fBx = fRx; fBy = fRy; + } + else + { + fAx = fRx; fAy = fRy; + } + // if last iteration brought too small advance, then do bisection next + // time, for safety + bHasToInterpolate = bHasToInterpolate && (std::abs(fRy) * 2.0 <= std::abs(fQy)); + ++nCount; + } + return fRx; +} + +// General functions + +void ScInterpreter::ScNoName() +{ + PushError(FormulaError::NoName); +} + +void ScInterpreter::ScBadName() +{ + short nParamCount = GetByte(); + while (nParamCount-- > 0) + { + PopError(); + } + PushError( FormulaError::NoName); +} + +double ScInterpreter::phi(double x) +{ + return 0.39894228040143268 * exp(-(x * x) / 2.0); +} + +double ScInterpreter::integralPhi(double x) +{ // Using gauss(x)+0.5 has severe cancellation errors for x<-4 + return 0.5 * std::erfc(-x * 0.7071067811865475); // * 1/sqrt(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 +/** you must ensure fZ>0 + Uses a variant of the Lanczos sum with a rational function. */ +static double lcl_getLanczosSum(double fZ) +{ + static const double fNum[13] ={ + 23531376880.41075968857200767445163675473, + 42919803642.64909876895789904700198885093, + 35711959237.35566804944018545154716670596, + 17921034426.03720969991975575445893111267, + 6039542586.35202800506429164430729792107, + 1439720407.311721673663223072794912393972, + 248874557.8620541565114603864132294232163, + 31426415.58540019438061423162831820536287, + 2876370.628935372441225409051620849613599, + 186056.2653952234950402949897160456992822, + 8071.672002365816210638002902272250613822, + 210.8242777515793458725097339207133627117, + 2.506628274631000270164908177133837338626 + }; + static const double fDenom[13] = { + 0, + 39916800, + 120543840, + 150917976, + 105258076, + 45995730, + 13339535, + 2637558, + 357423, + 32670, + 1925, + 66, + 1 + }; + // Horner scheme + double fSumNum; + double fSumDenom; + int nI; + if (fZ<=1.0) + { + fSumNum = fNum[12]; + fSumDenom = fDenom[12]; + for (nI = 11; nI >= 0; --nI) + { + fSumNum *= fZ; + fSumNum += fNum[nI]; + fSumDenom *= fZ; + fSumDenom += fDenom[nI]; + } + } + else + // Cancel down with fZ^12; Horner scheme with reverse coefficients + { + double fZInv = 1/fZ; + fSumNum = fNum[0]; + fSumDenom = fDenom[0]; + for (nI = 1; nI <=12; ++nI) + { + fSumNum *= fZInv; + fSumNum += fNum[nI]; + fSumDenom *= fZInv; + fSumDenom += fDenom[nI]; + } + } + return fSumNum/fSumDenom; +} + +// The algorithm is based on tgamma in gamma.hpp +// in math library from http://www.boost.org +/** You must ensure fZ>0; fZ>171.624376956302 will overflow. */ +static double lcl_GetGammaHelper(double fZ) +{ + double fGamma = lcl_getLanczosSum(fZ); + const double fg = 6.024680040776729583740234375; + double fZgHelp = fZ + fg - 0.5; + // avoid intermediate overflow + double fHalfpower = pow( fZgHelp, fZ / 2 - 0.25); + fGamma *= fHalfpower; + fGamma /= exp(fZgHelp); + fGamma *= fHalfpower; + if (fZ <= 20.0 && fZ == ::rtl::math::approxFloor(fZ)) + fGamma = ::rtl::math::round(fGamma); + return fGamma; +} + +// The algorithm is based on tgamma in gamma.hpp +// in math library from http://www.boost.org +/** You must ensure fZ>0 */ +static double lcl_GetLogGammaHelper(double fZ) +{ + const double fg = 6.024680040776729583740234375; + double fZgHelp = fZ + fg - 0.5; + return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp; +} + +/** You must ensure non integer arguments for fZ<1 */ +double ScInterpreter::GetGamma(double fZ) +{ + const double fLogPi = log(M_PI); + const double fLogDblMax = log( ::std::numeric_limits<double>::max()); + + if (fZ > fMaxGammaArgument) + { + SetError(FormulaError::IllegalFPOperation); + return HUGE_VAL; + } + + if (fZ >= 1.0) + return lcl_GetGammaHelper(fZ); + + if (fZ >= 0.5) // shift to x>=1 using Gamma(x)=Gamma(x+1)/x + return lcl_GetGammaHelper(fZ+1) / fZ; + + if (fZ >= -0.5) // shift to x>=1, might overflow + { + double fLogTest = lcl_GetLogGammaHelper(fZ+2) - rtl::math::log1p(fZ) - log( std::abs(fZ)); + if (fLogTest >= fLogDblMax) + { + SetError( FormulaError::IllegalFPOperation); + return HUGE_VAL; + } + return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ; + } + // fZ<-0.5 + // Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)*sin(pi*x) ) + double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( std::abs( ::rtl::math::sin( M_PI*fZ))); + if (fLogDivisor - fLogPi >= fLogDblMax) // underflow + return 0.0; + + if (fLogDivisor<0.0) + if (fLogPi - fLogDivisor > fLogDblMax) // overflow + { + SetError(FormulaError::IllegalFPOperation); + return HUGE_VAL; + } + + return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( M_PI*fZ) < 0.0) ? -1.0 : 1.0); +} + +/** You must ensure fZ>0 */ +double ScInterpreter::GetLogGamma(double fZ) +{ + if (fZ >= fMaxGammaArgument) + return lcl_GetLogGammaHelper(fZ); + if (fZ >= 1.0) + return log(lcl_GetGammaHelper(fZ)); + if (fZ >= 0.5) + return log( lcl_GetGammaHelper(fZ+1) / fZ); + return lcl_GetLogGammaHelper(fZ+2) - rtl::math::log1p(fZ) - log(fZ); +} + +double ScInterpreter::GetFDist(double x, double fF1, double fF2) +{ + double arg = fF2/(fF2+fF1*x); + double alpha = fF2/2.0; + double beta = fF1/2.0; + return GetBetaDist(arg, alpha, beta); +} + +double ScInterpreter::GetTDist( double T, double fDF, int nType ) +{ + switch ( nType ) + { + case 1 : // 1-tailed T-distribution + return 0.5 * GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5 ); + case 2 : // 2-tailed T-distribution + return GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5); + case 3 : // left-tailed T-distribution (probability density function) + return pow( 1 + ( T * T / fDF ), -( fDF + 1 ) / 2 ) / ( sqrt( fDF ) * GetBeta( 0.5, fDF / 2.0 ) ); + case 4 : // left-tailed T-distribution (cumulative distribution function) + double X = fDF / ( T * T + fDF ); + double R = 0.5 * GetBetaDist( X, 0.5 * fDF, 0.5 ); + return ( T < 0 ? R : 1 - R ); + } + SetError( FormulaError::IllegalArgument ); + return HUGE_VAL; +} + +// for LEGACY.CHIDIST, returns right tail, fDF=degrees of freedom +/** You must ensure fDF>0.0 */ +double ScInterpreter::GetChiDist(double fX, double fDF) +{ + if (fX <= 0.0) + return 1.0; // see ODFF + else + return GetUpRegIGamma( fDF/2.0, fX/2.0); +} + +// ready for ODF 1.2 +// for ODF CHISQDIST; cumulative distribution function, fDF=degrees of freedom +// returns left tail +/** You must ensure fDF>0.0 */ +double ScInterpreter::GetChiSqDistCDF(double fX, double fDF) +{ + if (fX <= 0.0) + return 0.0; // see ODFF + else + return GetLowRegIGamma( fDF/2.0, fX/2.0); +} + +double ScInterpreter::GetChiSqDistPDF(double fX, double fDF) +{ + // you must ensure fDF is positive integer + double fValue; + if (fX <= 0.0) + return 0.0; // see ODFF + if (fDF*fX > 1391000.0) + { + // intermediate invalid values, use log + fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0) - GetLogGamma(0.5*fDF)); + } + else // fDF is small in most cases, we can iterate + { + double fCount; + if (fmod(fDF,2.0)<0.5) + { + // even + fValue = 0.5; + fCount = 2.0; + } + else + { + fValue = 1/sqrt(fX*2*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 * ::rtl::math::log1p(fTempA) + -fB * ::rtl::math::log1p(fTempB)-fgm); + fResult *= fLanczos; + return fResult; +} + +// Same as GetBeta but with logarithm +double ScInterpreter::GetLogBeta(double fAlpha, double fBeta) +{ + double fA; + double fB; + if (fAlpha > fBeta) + { + fA = fAlpha; fB = fBeta; + } + else + { + fA = fBeta; fB = fAlpha; + } + const double fg = 6.024680040776729583740234375; //see GetGamma + double fgm = fg - 0.5; + double fLanczos = lcl_getLanczosSum(fA); + fLanczos /= lcl_getLanczosSum(fA+fB); + fLanczos *= lcl_getLanczosSum(fB); + double fLogLanczos = log(fLanczos); + double fABgm = fA+fB+fgm; + fLogLanczos += 0.5*(log(fABgm)-log(fA+fgm)-log(fB+fgm)); + double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm)) + double fTempB = fA/(fB+fgm); + double fResult = -fA * ::rtl::math::log1p(fTempA) + -fB * ::rtl::math::log1p(fTempB)-fgm; + fResult += fLogLanczos; + return fResult; +} + +// beta distribution probability density function +double ScInterpreter::GetBetaDistPDF(double fX, double fA, double fB) +{ + // special cases + if (fA == 1.0) // result b*(1-x)^(b-1) + { + if (fB == 1.0) + return 1.0; + if (fB == 2.0) + return -2.0*fX + 2.0; + if (fX == 1.0 && fB < 1.0) + { + SetError(FormulaError::IllegalArgument); + return HUGE_VAL; + } + if (fX <= 0.01) + return fB + fB * std::expm1((fB-1.0) * ::rtl::math::log1p(-fX)); + else + return fB * pow(0.5-fX+0.5,fB-1.0); + } + if (fB == 1.0) // result a*x^(a-1) + { + if (fA == 2.0) + return fA * fX; + if (fX == 0.0 && fA < 1.0) + { + SetError(FormulaError::IllegalArgument); + return HUGE_VAL; + } + return fA * pow(fX,fA-1); + } + if (fX <= 0.0) + { + if (fA < 1.0 && fX == 0.0) + { + SetError(FormulaError::IllegalArgument); + return HUGE_VAL; + } + else + return 0.0; + } + if (fX >= 1.0) + { + if (fB < 1.0 && fX == 1.0) + { + SetError(FormulaError::IllegalArgument); + return HUGE_VAL; + } + else + return 0.0; + } + + // normal cases; result x^(a-1)*(1-x)^(b-1)/Beta(a,b) + const double fLogDblMax = log( ::std::numeric_limits<double>::max()); + const double fLogDblMin = log( ::std::numeric_limits<double>::min()); + double fLogY = (fX < 0.1) ? ::rtl::math::log1p(-fX) : log(0.5-fX+0.5); + double fLogX = log(fX); + double fAm1LogX = (fA-1.0) * fLogX; + double fBm1LogY = (fB-1.0) * fLogY; + double fLogBeta = GetLogBeta(fA,fB); + // check whether parts over- or underflow + if ( fAm1LogX < fLogDblMax && fAm1LogX > fLogDblMin + && fBm1LogY < fLogDblMax && fBm1LogY > fLogDblMin + && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin + && fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > fLogDblMin) + return pow(fX,fA-1.0) * pow(0.5-fX+0.5,fB-1.0) / GetBeta(fA,fB); + else // need logarithm; + // might overflow as a whole, but seldom, not worth to pre-detect it + return exp( fAm1LogX + fBm1LogY - fLogBeta); +} + +/* + x^a * (1-x)^b + I_x(a,b) = ---------------- * result of ContFrac + a * Beta(a,b) +*/ +static double lcl_GetBetaHelperContFrac(double fX, double fA, double fB) +{ // like old version + double a1, b1, a2, b2, fnorm, cfnew, cf; + a1 = 1.0; b1 = 1.0; + b2 = 1.0 - (fA+fB)/(fA+1.0)*fX; + if (b2 == 0.0) + { + a2 = 0.0; + fnorm = 1.0; + cf = 1.0; + } + else + { + a2 = 1.0; + fnorm = 1.0/b2; + cf = a2*fnorm; + } + cfnew = 1.0; + double rm = 1.0; + + const double fMaxIter = 50000.0; + // loop security, normal cases converge in less than 100 iterations. + // FIXME: You will get so much iterations for fX near mean, + // I do not know a better algorithm. + bool bfinished = false; + do + { + const double apl2m = fA + 2.0*rm; + const double d2m = rm*(fB-rm)*fX/((apl2m-1.0)*apl2m); + const double d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+1.0)); + a1 = (a2+d2m*a1)*fnorm; + b1 = (b2+d2m*b1)*fnorm; + a2 = a1 + d2m1*a2*fnorm; + b2 = b1 + d2m1*b2*fnorm; + if (b2 != 0.0) + { + fnorm = 1.0/b2; + cfnew = a2*fnorm; + bfinished = (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 * ::rtl::math::log1p(-fXin)); + //FIXME: need special algorithm for fX near fP for large fA,fB + double fResult; + // I use always continued fraction, power series are neither + // faster nor more accurate. + double fY = (0.5-fXin)+0.5; + double flnY = ::rtl::math::log1p(-fXin); + double fX = fXin; + double flnX = log(fXin); + double fA = fAlpha; + double fB = fBeta; + bool bReflect = fXin > fAlpha/(fAlpha+fBeta); + if (bReflect) + { + fA = fBeta; + fB = fAlpha; + fX = fY; + fY = fXin; + flnX = flnY; + flnY = log(fXin); + } + fResult = lcl_GetBetaHelperContFrac(fX,fA,fB); + fResult = fResult/fA; + double fP = fA/(fA+fB); + double fQ = fB/(fA+fB); + double fTemp; + if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97) //found experimental + fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY; + else + fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB)); + fResult *= fTemp; + if (bReflect) + fResult = 0.5 - fResult + 0.5; + if (fResult > 1.0) // ensure valid range + fResult = 1.0; + if (fResult < 0.0) + fResult = 0.0; + return fResult; +} + +void ScInterpreter::ScBetaDist() +{ + sal_uInt8 nParamCount = GetByte(); + if ( !MustHaveParamCount( nParamCount, 3, 6 ) ) // expanded, see #i91547# + return; + double fLowerBound, fUpperBound; + double alpha, beta, x; + bool bIsCumulative; + if (nParamCount == 6) + bIsCumulative = GetBool(); + else + bIsCumulative = true; + if (nParamCount >= 5) + fUpperBound = GetDouble(); + else + fUpperBound = 1.0; + if (nParamCount >= 4) + fLowerBound = GetDouble(); + else + fLowerBound = 0.0; + beta = GetDouble(); + alpha = GetDouble(); + x = GetDouble(); + double fScale = fUpperBound - fLowerBound; + if (fScale <= 0.0 || alpha <= 0.0 || beta <= 0.0) + { + PushIllegalArgument(); + return; + } + if (bIsCumulative) // cumulative distribution function + { + // special cases + if (x < fLowerBound) + { + PushDouble(0.0); return; //see spec + } + if (x > fUpperBound) + { + PushDouble(1.0); return; //see spec + } + // normal cases + x = (x-fLowerBound)/fScale; // convert to standard form + PushDouble(GetBetaDist(x, alpha, beta)); + return; + } + else // probability density function + { + if (x < fLowerBound || x > fUpperBound) + { + PushDouble(0.0); + return; + } + x = (x-fLowerBound)/fScale; + PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale); + return; + } +} + +/** + Microsoft version has parameters in different order + Also, upper and lowerbound are optional and have default values + and different constraints apply. + Basically, function is identical with ScInterpreter::ScBetaDist() +*/ +void ScInterpreter::ScBetaDist_MS() +{ + sal_uInt8 nParamCount = GetByte(); + if ( !MustHaveParamCount( nParamCount, 4, 6 ) ) + return; + double fLowerBound, fUpperBound; + double alpha, beta, x; + bool bIsCumulative; + if (nParamCount == 6) + fUpperBound = GetDouble(); + else + fUpperBound = 1.0; + if (nParamCount >= 5) + fLowerBound = GetDouble(); + else + fLowerBound = 0.0; + bIsCumulative = GetBool(); + beta = GetDouble(); + alpha = GetDouble(); + x = GetDouble(); + if (alpha <= 0.0 || beta <= 0.0 || x < fLowerBound || x > fUpperBound) + { + PushIllegalArgument(); + return; + } + double fScale = fUpperBound - fLowerBound; + if (bIsCumulative) // cumulative distribution function + { + x = (x-fLowerBound)/fScale; // convert to standard form + PushDouble(GetBetaDist(x, alpha, beta)); + return; + } + else // probability density function + { + x = (x-fLowerBound)/fScale; + PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale); + return; + } +} + +void ScInterpreter::ScPhi() +{ + PushDouble(phi(GetDouble())); +} + +void ScInterpreter::ScGauss() +{ + PushDouble(gauss(GetDouble())); +} + +void ScInterpreter::ScFisher() +{ + double fVal = GetDouble(); + if (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()); + } + } + } +} + +/** Local function used in the calculation of the hypergeometric distribution. + */ +static void lcl_PutFactorialElements( ::std::vector< double >& cn, double fLower, double fUpper, double fBase ) +{ + for ( double i = fLower; i <= fUpper; ++i ) + { + double fVal = fBase - i; + if ( fVal > 1.0 ) + cn.push_back( fVal ); + } +} + +/** Calculates a value of the hypergeometric distribution. + + @see #i47296# + + This function has an extra argument bCumulative, + which only calculates the non-cumulative distribution and + which is optional in Calc and mandatory with Excel's HYPGEOM.DIST() + + @see fdo#71722 + @see tdf#102948, make Calc function ODFF1.2-compliant + @see tdf#117041, implement note at bottom of ODFF1.2 par.6.18.37 + */ +void ScInterpreter::ScHypGeomDist( int nMinParamCount ) +{ + sal_uInt8 nParamCount = GetByte(); + if ( !MustHaveParamCount( nParamCount, nMinParamCount, 5 ) ) + return; + + bool bCumulative = ( nParamCount == 5 && GetBool() ); + double N = ::rtl::math::approxFloor(GetDouble()); + double M = ::rtl::math::approxFloor(GetDouble()); + double n = ::rtl::math::approxFloor(GetDouble()); + double x = ::rtl::math::approxFloor(GetDouble()); + + if ( (x < 0.0) || (n < x) || (N < n) || (N < M) || (M < 0.0) ) + { + PushIllegalArgument(); + return; + } + + KahanSum fVal = 0.0; + + for ( int i = ( bCumulative ? 0 : x ); i <= x && nGlobalError == FormulaError::NONE; i++ ) + { + if ( (n - i <= N - M) && (i <= M) ) + fVal += GetHypGeomDist( i, n, M, N ); + } + + PushDouble( fVal.get() ); +} + +/** Calculates a value of the hypergeometric distribution. + + The algorithm is designed to avoid unnecessary multiplications and division + by expanding all factorial elements (9 of them). It is done by excluding + those ranges that overlap in the numerator and the denominator. This allows + for a fast calculation for large values which would otherwise cause an overflow + in the intermediate values. + + @see #i47296# + */ +double ScInterpreter::GetHypGeomDist( double x, double n, double M, double N ) +{ + const size_t nMaxArraySize = 500000; // arbitrary max array size + + std::vector<double> cnNumer, cnDenom; + + size_t nEstContainerSize = static_cast<size_t>( x + ::std::min( n, M ) ); + size_t nMaxSize = ::std::min( cnNumer.max_size(), nMaxArraySize ); + if ( nEstContainerSize > nMaxSize ) + { + PushNoValue(); + return 0; + } + cnNumer.reserve( nEstContainerSize + 10 ); + cnDenom.reserve( nEstContainerSize + 10 ); + + // Trim coefficient C first + double fCNumVarUpper = N - n - M + x - 1.0; + double fCDenomVarLower = 1.0; + if ( N - n - M + x >= M - x + 1.0 ) + { + fCNumVarUpper = M - x - 1.0; + fCDenomVarLower = N - n - 2.0*(M - x) + 1.0; + } + + double fCNumLower = N - n - fCNumVarUpper; + double fCDenomUpper = N - n - M + x + 1.0 - fCDenomVarLower; + + double fDNumVarLower = n - M; + + if ( n >= M + 1.0 ) + { + if ( N - M < n + 1.0 ) + { + // Case 1 + + if ( N - n < n + 1.0 ) + { + // no overlap + lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n ); + lcl_PutFactorialElements( cnDenom, 0.0, N - n - 1.0, N ); + } + else + { + // overlap + OSL_ENSURE( fCNumLower < n + 1.0, "ScHypGeomDist: wrong assertion" ); + lcl_PutFactorialElements( cnNumer, N - 2.0*n, fCNumVarUpper, N - n ); + lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N ); + } + + OSL_ENSURE( fCDenomUpper <= N - M, "ScHypGeomDist: wrong assertion" ); + + if ( fCDenomUpper < n - x + 1.0 ) + // no overlap + lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 ); + else + { + // overlap + lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 ); + + fCDenomUpper = n - x; + fCDenomVarLower = N - M - 2.0*(n - x) + 1.0; + } + } + else + { + // Case 2 + + if ( n > M - 1.0 ) + { + // no overlap + lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n ); + lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N ); + } + else + { + lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n ); + lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N ); + } + + OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" ); + + if ( fCDenomUpper < n - x + 1.0 ) + // no overlap + lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - n + x, N - M + 1.0 ); + else + { + lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - fCDenomUpper, N - M + 1.0 ); + fCDenomUpper = n - x; + fCDenomVarLower = N - M - 2.0*(n - x) + 1.0; + } + } + + OSL_ENSURE( fCDenomUpper <= M, "ScHypGeomDist: wrong assertion" ); + } + else + { + if ( N - M < M + 1.0 ) + { + // Case 3 + + if ( N - n < M + 1.0 ) + { + // No overlap + lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n ); + lcl_PutFactorialElements( cnDenom, 0.0, N - M - 1.0, N ); + } + else + { + lcl_PutFactorialElements( cnNumer, N - n - M, fCNumVarUpper, N - n ); + lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N ); + } + + if ( n - x + 1.0 > fCDenomUpper ) + // No overlap + lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 ); + else + { + // Overlap + lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 ); + + fCDenomVarLower = N - M - 2.0*(n - x) + 1.0; + fCDenomUpper = n - x; + } + } + else + { + // Case 4 + + OSL_ENSURE( M >= n - x, "ScHypGeomDist: wrong assertion" ); + OSL_ENSURE( M - x <= N - M + 1.0, "ScHypGeomDist: wrong assertion" ); + + if ( N - n < N - M + 1.0 ) + { + // No overlap + lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n ); + lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N ); + } + else + { + // Overlap + OSL_ENSURE( fCNumLower <= N - M + 1.0, "ScHypGeomDist: wrong assertion" ); + lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n ); + lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N ); + } + + if ( n - x + 1.0 > fCDenomUpper ) + // No overlap + lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - n + x, N - M + 1.0 ); + else if ( M >= fCDenomUpper ) + { + lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - fCDenomUpper, N - M + 1.0 ); + + fCDenomUpper = n - x; + fCDenomVarLower = N - M - 2.0*(n - x) + 1.0; + } + else + { + OSL_ENSURE( M <= fCDenomUpper, "ScHypGeomDist: wrong assertion" ); + lcl_PutFactorialElements( cnDenom, fCDenomVarLower, N - n - 2.0*M + x, + N - n - M + x + 1.0 ); + + fCDenomUpper = n - x; + fCDenomVarLower = N - M - 2.0*(n - x) + 1.0; + } + } + + OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" ); + + fDNumVarLower = 0.0; + } + + double nDNumVarUpper = fCDenomUpper < x + 1.0 ? n - x - 1.0 : n - fCDenomUpper - 1.0; + double nDDenomVarLower = fCDenomUpper < x + 1.0 ? fCDenomVarLower : N - n - M + 1.0; + lcl_PutFactorialElements( cnNumer, fDNumVarLower, nDNumVarUpper, n ); + lcl_PutFactorialElements( cnDenom, nDDenomVarLower, N - n - M + x, N - n - M + x + 1.0 ); + + ::std::sort( cnNumer.begin(), cnNumer.end() ); + ::std::sort( cnDenom.begin(), cnDenom.end() ); + auto it1 = cnNumer.rbegin(), it1End = cnNumer.rend(); + auto it2 = cnDenom.rbegin(), it2End = cnDenom.rend(); + + double fFactor = 1.0; + for ( ; it1 != it1End || it2 != it2End; ) + { + double fEnum = 1.0, fDenom = 1.0; + if ( it1 != it1End ) + fEnum = *it1++; + if ( it2 != it2End ) + fDenom = *it2++; + fFactor *= fEnum / fDenom; + } + + return fFactor; +} + +void ScInterpreter::ScGammaDist( bool bODFF ) +{ + sal_uInt8 nMinParamCount = ( bODFF ? 3 : 4 ); + sal_uInt8 nParamCount = GetByte(); + if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) ) + return; + bool bCumulative; + if (nParamCount == 4) + bCumulative = GetBool(); + else + bCumulative = true; + double fBeta = GetDouble(); // scale + double fAlpha = GetDouble(); // shape + double fX = GetDouble(); // x + if ((!bODFF && fX < 0) || fAlpha <= 0.0 || fBeta <= 0.0) + PushIllegalArgument(); + else + { + if (bCumulative) // distribution + PushDouble( GetGammaDist( fX, fAlpha, fBeta)); + else // density + PushDouble( GetGammaDistPDF( fX, fAlpha, fBeta)); + } +} + +void ScInterpreter::ScNormInv() +{ + if ( MustHaveParamCount( GetByte(), 3 ) ) + { + double sigma = GetDouble(); + double mue = GetDouble(); + double x = GetDouble(); + if (sigma <= 0.0 || x < 0.0 || x > 1.0) + PushIllegalArgument(); + else if (x == 0.0 || x == 1.0) + PushNoValue(); + else + PushDouble(gaussinv(x)*sigma + mue); + } +} + +void ScInterpreter::ScSNormInv() +{ + double x = GetDouble(); + if (x < 0.0 || x > 1.0) + PushIllegalArgument(); + else if (x == 0.0 || x == 1.0) + PushNoValue(); + else + PushDouble(gaussinv(x)); +} + +void ScInterpreter::ScLogNormInv() +{ + sal_uInt8 nParamCount = GetByte(); + if ( MustHaveParamCount( nParamCount, 1, 3 ) ) + { + double fSigma = ( nParamCount == 3 ? GetDouble() : 1.0 ); // Stddev + double fMue = ( nParamCount >= 2 ? GetDouble() : 0.0 ); // Mean + double fP = GetDouble(); // p + if ( fSigma <= 0.0 || fP <= 0.0 || fP >= 1.0 ) + PushIllegalArgument(); + else + PushDouble( exp( fMue + fSigma * gaussinv( fP ) ) ); + } +} + +class ScGammaDistFunction : public ScDistFunc +{ + ScInterpreter& rInt; + double fp, fAlpha, fBeta; + +public: + ScGammaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) : + rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {} + + virtual ~ScGammaDistFunction() {} + + double GetValue( double x ) const override { return fp - rInt.GetGammaDist(x, fAlpha, fBeta); } +}; + +void ScInterpreter::ScGammaInv() +{ + if ( !MustHaveParamCount( GetByte(), 3 ) ) + return; + double fBeta = GetDouble(); + double fAlpha = GetDouble(); + 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, mrDoc, 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, mrDoc, 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, mrDoc, 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, mrDoc, 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 = GetTopNumberArray(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), + [](double f) { + f = rtl::math::approxFloor(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::GetTopNumberArray( 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, mrDoc, aRange, mnSubTotalFlags); + if (aValIter.GetFirst(fCellVal, nErr)) + { + do + aArray.push_back(fCellVal); + while (aValIter.GetNext(fCellVal, nErr) && nErr == FormulaError::NONE); + } + 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, mrDoc, 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 ] ) + PushNoValue(); + 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 ( !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, mrDoc, 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, mrDoc, 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) + { + double fTemp = fLo; + fLo = fUp; + fUp = fTemp; + } + ScMatrixRef pMatP = GetMatrix(); + ScMatrixRef pMatW = GetMatrix(); + if (!pMatP || !pMatW) + PushIllegalParameter(); + else + { + SCSIZE nC1, nC2; + SCSIZE nR1, nR2; + pMatP->GetDimensions(nC1, nR1); + pMatW->GetDimensions(nC2, nR2); + if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 || + nC2 == 0 || nR2 == 0) + PushNA(); + else + { + 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 = sqrt(fR*fR + fI*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: */ |