diff options
Diffstat (limited to 'sc/source/core/tool/interpr3.cxx')
| -rw-r--r-- | sc/source/core/tool/interpr3.cxx | 4157 |
1 files changed, 0 insertions, 4157 deletions
diff --git a/sc/source/core/tool/interpr3.cxx b/sc/source/core/tool/interpr3.cxx deleted file mode 100644 index 62d670d65..000000000 --- a/sc/source/core/tool/interpr3.cxx +++ /dev/null @@ -1,4157 +0,0 @@ -/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ -/************************************************************************* - * - * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. - * - * Copyright 2000, 2010 Oracle and/or its affiliates. - * - * OpenOffice.org - a multi-platform office productivity suite - * - * This file is part of OpenOffice.org. - * - * OpenOffice.org is free software: you can redistribute it and/or modify - * it under the terms of the GNU Lesser General Public License version 3 - * only, as published by the Free Software Foundation. - * - * OpenOffice.org is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU Lesser General Public License version 3 for more details - * (a copy is included in the LICENSE file that accompanied this code). - * - * You should have received a copy of the GNU Lesser General Public License - * version 3 along with OpenOffice.org. If not, see - * <http://www.openoffice.org/license.html> - * for a copy of the LGPLv3 License. - * - ************************************************************************/ - -// MARKER(update_precomp.py): autogen include statement, do not remove -#include "precompiled_sc.hxx" - -// INCLUDE --------------------------------------------------------------- - -#include <tools/solar.h> -#include <stdlib.h> -#include <string.h> -#include <rtl/logfile.hxx> - -#include "interpre.hxx" -#include "global.hxx" -#include "compiler.hxx" -#include "cell.hxx" -#include "document.hxx" -#include "dociter.hxx" -#include "scmatrix.hxx" -#include "globstr.hrc" - -#include <math.h> -#include <vector> -#include <algorithm> - -using ::std::vector; -using namespace formula; - -// STATIC DATA ----------------------------------------------------------- - -#define SCdEpsilon 1.0E-7 -#define SC_MAX_ITERATION_COUNT 20 -#define MAX_ANZ_DOUBLE_FOR_SORT 100000 - -const double ScInterpreter::fMaxGammaArgument = 171.624376956302; // found experimental -const double fMachEps = ::std::numeric_limits<double>::epsilon(); - -//----------------------------------------------------------------------------- - -class ScDistFunc -{ -public: - virtual double GetValue(double x) const = 0; -}; - -// 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 */ -inline bool lcl_HasChangeOfSign( double u, double w ) -{ - return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0); -} - -double lcl_IterateInverse( const ScDistFunc& rFunction, double fAx, double fBx, bool& rConvError ) -{ - rConvError = false; - const double fYEps = 1.0E-307; - const double fXEps = ::std::numeric_limits<double>::epsilon(); - - OSL_ENSURE(fAx<fBx, "IterateInverse: wrong interval"); - - // find enclosing interval - - double fAy = rFunction.GetValue(fAx); - double fBy = rFunction.GetValue(fBx); - double fTemp; - unsigned short nCount; - for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy); nCount++) - { - if (fabs(fAy) <= fabs(fBy)) - { - fTemp = fAx; - fAx += 2.0 * (fAx - fBx); - if (fAx < 0.0) - fAx = 0.0; - fBx = fTemp; - fBy = fAy; - fAy = rFunction.GetValue(fAx); - } - else - { - fTemp = fBx; - fBx += 2.0 * (fBx - fAx); - fAx = fTemp; - fAy = fBy; - fBy = rFunction.GetValue(fBx); - } - } - - if (fAy == 0.0) - return fAx; - if (fBy == 0.0) - return fBx; - if (!lcl_HasChangeOfSign( fAy, fBy)) - { - rConvError = true; - return 0.0; - } - // inverse quadric interpolation with additional brackets - // set three points - double fPx = fAx; - double fPy = fAy; - double fQx = fBx; - double fQy = fBy; - double fRx = fAx; - double fRy = fAy; - double fSx = 0.5 * (fAx + fBx); // potential next point - bool bHasToInterpolate = true; - nCount = 0; - while ( nCount < 500 && fabs(fRy) > fYEps && - (fBx-fAx) > ::std::max( fabs(fAx), fabs(fBx)) * fXEps ) - { - if (bHasToInterpolate) - { - if (fPy!=fQy && fQy!=fRy && fRy!=fPy) - { - fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy) - + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy) - + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy); - bHasToInterpolate = (fAx < fSx) && (fSx < fBx); // inside the brackets? - } - else - bHasToInterpolate = false; - } - if(!bHasToInterpolate) - { - fSx = 0.5 * (fAx + fBx); - // reset points - fPx = fAx; fPy = fAy; - 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 interration brought to small advance, then do bisection next - // time, for safety - bHasToInterpolate = bHasToInterpolate && (fabs(fRy) * 2.0 <= fabs(fQy)); - ++nCount; - } - return fRx; -} - -//----------------------------------------------------------------------------- -// Allgemeine Funktionen -//----------------------------------------------------------------------------- - -void ScInterpreter::ScNoName() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScNoName" ); - PushError(errNoName); -} - -void ScInterpreter::ScBadName() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScBadName" ); - short nParamCount = GetByte(); - while (nParamCount-- > 0) - { - PopError(); - } - PushError( errNoName); -} - -double ScInterpreter::phi(double x) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::phi" ); - return 0.39894228040143268 * exp(-(x * x) / 2.0); -} - -double ScInterpreter::integralPhi(double x) -{ // Using gauss(x)+0.5 has severe cancellation errors for x<-4 - return 0.5 * ::rtl::math::erfc(-x * 0.7071067811865475); // * 1/sqrt(2) -} - -double ScInterpreter::taylor(double* pPolynom, sal_uInt16 nMax, double x) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::taylor" ); - double nVal = pPolynom[nMax]; - for (short i = nMax-1; i >= 0; i--) - { - nVal = pPolynom[i] + (nVal * x); - } - return nVal; -} - -double ScInterpreter::gauss(double x) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::gauss" ); - 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 }; - 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 }; - 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 }; - double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 }; - - double xAbs = fabs(x); - sal_uInt16 xShort = (sal_uInt16)::rtl::math::approxFloor(xAbs); - double nVal = 0.0; - if (xShort == 0) - nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs; - else if ((xShort >= 1) && (xShort <= 2)) - nVal = taylor(t2, 23, (xAbs - 2.0)); - else if ((xShort >= 3) && (xShort <= 4)) - nVal = taylor(t4, 20, (xAbs - 4.0)); - else - 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) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::gaussinv" ); - double q,t,z; - - q=x-0.5; - - if(fabs(q)<=.425) - { - t=0.180625-q*q; - - z= - q* - ( - ( - ( - ( - ( - ( - ( - t*2509.0809287301226727+33430.575583588128105 - ) - *t+67265.770927008700853 - ) - *t+45921.953931549871457 - ) - *t+13731.693765509461125 - ) - *t+1971.5909503065514427 - ) - *t+133.14166789178437745 - ) - *t+3.387132872796366608 - ) - / - ( - ( - ( - ( - ( - ( - ( - t*5226.495278852854561+28729.085735721942674 - ) - *t+39307.89580009271061 - ) - *t+21213.794301586595867 - ) - *t+5394.1960214247511077 - ) - *t+687.1870074920579083 - ) - *t+42.313330701600911252 - ) - *t+1.0 - ); - - } - else - { - if(q>0) t=1-x; - else t=x; - - t=sqrt(-log(t)); - - if(t<=5.0) - { - t+=-1.6; - - z= - ( - ( - ( - ( - ( - ( - ( - t*7.7454501427834140764e-4+0.0227238449892691845833 - ) - *t+0.24178072517745061177 - ) - *t+1.27045825245236838258 - ) - *t+3.64784832476320460504 - ) - *t+5.7694972214606914055 - ) - *t+4.6303378461565452959 - ) - *t+1.42343711074968357734 - ) - / - ( - ( - ( - ( - ( - ( - ( - t*1.05075007164441684324e-9+5.475938084995344946e-4 - ) - *t+0.0151986665636164571966 - ) - *t+0.14810397642748007459 - ) - *t+0.68976733498510000455 - ) - *t+1.6763848301838038494 - ) - *t+2.05319162663775882187 - ) - *t+1.0 - ); - - } - else - { - t+=-5.0; - - z= - ( - ( - ( - ( - ( - ( - ( - t*2.01033439929228813265e-7+2.71155556874348757815e-5 - ) - *t+0.0012426609473880784386 - ) - *t+0.026532189526576123093 - ) - *t+0.29656057182850489123 - ) - *t+1.7848265399172913358 - ) - *t+5.4637849111641143699 - ) - *t+6.6579046435011037772 - ) - / - ( - ( - ( - ( - ( - ( - ( - t*2.04426310338993978564e-15+1.4215117583164458887e-7 - ) - *t+1.8463183175100546818e-5 - ) - *t+7.868691311456132591e-4 - ) - *t+0.0148753612908506148525 - ) - *t+0.13692988092273580531 - ) - *t+0.59983220655588793769 - ) - *t+1.0 - ); - - } - - if(q<0.0) z=-z; - } - - return z; -} - -double ScInterpreter::Fakultaet(double x) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::Fakultaet" ); - 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(errNoValue); - return x; -} - -double ScInterpreter::BinomKoeff(double n, double k) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::BinomKoeff" ); - 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. */ -double lcl_getLanczosSum(double fZ) -{ - 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 - }; - 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. */ -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 */ -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) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGamma" ); - const double fLogPi = log(F_PI); - const double fLogDblMax = log( ::std::numeric_limits<double>::max()); - - if (fZ > fMaxGammaArgument) - { - SetError(errIllegalFPOperation); - 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) - log(fZ+1) - log( fabs(fZ)); - if (fLogTest >= fLogDblMax) - { - SetError( errIllegalFPOperation); - return HUGE_VAL; - } - return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ; - } - // fZ<-0.5 - // Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)*sin(pi*x) ) - double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( fabs( ::rtl::math::sin( F_PI*fZ))); - if (fLogDivisor - fLogPi >= fLogDblMax) // underflow - return 0.0; - - if (fLogDivisor<0.0) - if (fLogPi - fLogDivisor > fLogDblMax) // overflow - { - SetError(errIllegalFPOperation); - return HUGE_VAL; - } - - return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( F_PI*fZ) < 0.0) ? -1.0 : 1.0); -} - -/** You must ensure fZ>0 */ -double ScInterpreter::GetLogGamma(double fZ) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetLogGamma" ); - 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) - log(fZ+1) - log(fZ); -} - -double ScInterpreter::GetFDist(double x, double fF1, double fF2) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetFDist" ); - 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) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetTDist" ); - return 0.5 * GetBetaDist(fDF/(fDF+T*T), fDF/2.0, 0.5); -} - -// for LEGACY.CHIDIST, returns right tail, fDF=degrees of freedom -/** You must ensure fDF>0.0 */ -double ScInterpreter::GetChiDist(double fX, double fDF) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetChiDist" ); - 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) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetChiSqDistCDF" ); - if (fX <= 0.0) - return 0.0; // see ODFF - else - return GetLowRegIGamma( fDF/2.0, fX/2.0); -} - -double ScInterpreter::GetChiSqDistPDF(double fX, double fDF) -{ - // you must ensure fDF is positive integer - double fValue; - if (fX <= 0.0) - return 0.0; // see ODFF - if (fDF*fX > 1391000.0) - { - // intermediate invalid values, use log - fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0) - GetLogGamma(0.5*fDF)); - } - else // fDF is small in most cases, we can iterate - { - double fCount; - if (fmod(fDF,2.0)<0.5) - { - // even - fValue = 0.5; - fCount = 2.0; - } - else - { - fValue = 1/sqrt(fX*2*F_PI); - fCount = 1.0; - } - while ( fCount < fDF) - { - fValue *= (fX / fCount); - fCount += 2.0; - } - if (fX>=1425.0) // underflow in e^(-x/2) - fValue = exp(log(fValue)-fX/2); - else - fValue *= exp(-fX/2); - } - return fValue; -} - -void ScInterpreter::ScChiSqDist() -{ - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 2, 3 ) ) - return; - bool bCumulative; - if (nParamCount == 3) - bCumulative = GetBool(); - else - bCumulative = true; - double fDF = ::rtl::math::approxFloor(GetDouble()); - if (fDF < 1.0) - PushIllegalArgument(); - else - { - double fX = GetDouble(); - if (bCumulative) - PushDouble(GetChiSqDistCDF(fX,fDF)); - else - PushDouble(GetChiSqDistPDF(fX,fDF)); - } -} - -void ScInterpreter::ScGamma() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGamma" ); - double x = GetDouble(); - if (x <= 0.0 && x == ::rtl::math::approxFloor(x)) - PushIllegalArgument(); - else - { - double fResult = GetGamma(x); - if (nGlobalError) - { - PushError( nGlobalError); - return; - } - PushDouble(fResult); - } -} - -void ScInterpreter::ScLogGamma() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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(errIllegalArgument); - return HUGE_VAL; - } - if (fX <= 0.01) - return fB + fB * ::rtl::math::expm1((fB-1.0) * ::rtl::math::log1p(-fX)); - else - return fB * pow(0.5-fX+0.5,fB-1.0); - } - if (fB == 1.0) // result a*x^(a-1) - { - if (fA == 2.0) - return fA * fX; - if (fX == 0.0 && fA < 1.0) - { - SetError(errIllegalArgument); - return HUGE_VAL; - } - return fA * pow(fX,fA-1); - } - if (fX <= 0.0) - { - if (fA < 1.0 && fX == 0.0) - { - SetError(errIllegalArgument); - return HUGE_VAL; - } - else - return 0.0; - } - if (fX >= 1.0) - { - if (fB < 1.0 && fX == 1.0) - { - SetError(errIllegalArgument); - 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 fAm1 = fA-1.0; - double fBm1 = fB-1.0; - double fLogBeta = GetLogBeta(fA,fB); - // check whether parts over- or underflow - if ( fAm1 * fLogX < fLogDblMax && fAm1 * fLogX > fLogDblMin - && fBm1 * fLogY < fLogDblMax && fBm1* fLogY > fLogDblMin - && fLogBeta < fLogDblMax && fLogBeta > 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((fA-1.0)*fLogX + (fB-1.0)* fLogY - fLogBeta); -} - -/* - x^a * (1-x)^b - I_x(a,b) = ---------------- * result of ContFrac - a * Beta(a,b) -*/ -double lcl_GetBetaHelperContFrac(double fX, double fA, double fB) -{ // like old version - double a1, b1, a2, b2, fnorm, apl2m, d2m, d2m1, 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 iteratons for fX near mean, - // I do not know a better algorithm. - bool bfinished = false; - do - { - apl2m = fA + 2.0*rm; - d2m = rm*(fB-rm)*fX/((apl2m-1.0)*apl2m); - d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+1.0)); - a1 = (a2+d2m*a1)*fnorm; - b1 = (b2+d2m*b1)*fnorm; - a2 = a1 + d2m1*a2*fnorm; - b2 = b1 + d2m1*b2*fnorm; - if (b2 != 0.0) - { - fnorm = 1.0/b2; - cfnew = a2*fnorm; - bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps); - } - cf = cfnew; - rm += 1.0; - } - while (rm < fMaxIter && !bfinished); - return cf; -} - -// cumulative distribution function, normalized -double ScInterpreter::GetBetaDist(double fXin, double fAlpha, double fBeta) -{ -// special cases - if (fXin <= 0.0) // values are valid, see spec - return 0.0; - if (fXin >= 1.0) // values are valid, see spec - return 1.0; - if (fBeta == 1.0) - return pow(fXin, fAlpha); - if (fAlpha == 1.0) - // 1.0 - pow(1.0-fX,fBeta) is not accurate enough - return -::rtl::math::expm1(fBeta * ::rtl::math::log1p(-fXin)); - //FIXME: need special algorithm for fX near fP for large fA,fB - double fResult; - // I use always continued fraction, power series are neither - // faster nor more accurate. - double fY = (0.5-fXin)+0.5; - double flnY = ::rtl::math::log1p(-fXin); - double fX = fXin; - double flnX = log(fXin); - double fA = fAlpha; - double fB = fBeta; - bool bReflect = fXin > fAlpha/(fAlpha+fBeta); - if (bReflect) - { - fA = fBeta; - fB = fAlpha; - fX = fY; - fY = fXin; - flnX = flnY; - flnY = log(fXin); - } - fResult = lcl_GetBetaHelperContFrac(fX,fA,fB); - fResult = fResult/fA; - double fP = fA/(fA+fB); - double fQ = fB/(fA+fB); - double fTemp; - if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97) //found experimental - fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY; - else - fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB)); - fResult *= fTemp; - if (bReflect) - fResult = 0.5 - fResult + 0.5; - if (fResult > 1.0) // ensure valid range - fResult = 1.0; - if (fResult < 0.0) - fResult = 0.0; - return fResult; -} - - void ScInterpreter::ScBetaDist() - { - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 3, 6 ) ) // expanded, see #i91547# - return; - double fLowerBound, fUpperBound; - double alpha, beta, x; - bool bIsCumulative; - if (nParamCount == 6) - bIsCumulative = GetBool(); - else - bIsCumulative = true; - if (nParamCount >= 5) - fUpperBound = GetDouble(); - else - fUpperBound = 1.0; - if (nParamCount >= 4) - fLowerBound = GetDouble(); - else - fLowerBound = 0.0; - beta = GetDouble(); - alpha = GetDouble(); - x = GetDouble(); - double fScale = fUpperBound - fLowerBound; - if (fScale <= 0.0 || alpha <= 0.0 || beta <= 0.0) - { - PushIllegalArgument(); - return; - } - if (bIsCumulative) // cumulative distribution function - { - // special cases - if (x < fLowerBound) - { - PushDouble(0.0); return; //see spec - } - if (x > fUpperBound) - { - PushDouble(1.0); return; //see spec - } - // normal cases - x = (x-fLowerBound)/fScale; // convert to standard form - PushDouble(GetBetaDist(x, alpha, beta)); - return; - } - else // probability density function - { - if (x < fLowerBound || x > fUpperBound) - { - PushDouble(0.0); - return; - } - x = (x-fLowerBound)/fScale; - PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale); - return; - } -} - -void ScInterpreter::ScPhi() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScPhi" ); - PushDouble(phi(GetDouble())); -} - -void ScInterpreter::ScGauss() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGauss" ); - PushDouble(gauss(GetDouble())); -} - -void ScInterpreter::ScFisher() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScFisher" ); - double fVal = GetDouble(); - if (fabs(fVal) >= 1.0) - PushIllegalArgument(); - else - PushDouble( ::rtl::math::atanh( fVal)); -} - -void ScInterpreter::ScFisherInv() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScFisherInv" ); - PushDouble( tanh( GetDouble())); -} - -void ScInterpreter::ScFact() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScFact" ); - double nVal = GetDouble(); - if (nVal < 0.0) - PushIllegalArgument(); - else - PushDouble(Fakultaet(nVal)); -} - -void ScInterpreter::ScKombin() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScKombin" ); - 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::ScKombin2() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScKombin2" ); - 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::ScVariationen() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScVariationen" ); - if ( MustHaveParamCount( GetByte(), 2 ) ) - { - double k = ::rtl::math::approxFloor(GetDouble()); - double n = ::rtl::math::approxFloor(GetDouble()); - if (n < 0.0 || k < 0.0 || k > n) - PushIllegalArgument(); - else if (k == 0.0) - PushInt(1); // (n! / (n - 0)!) == 1 - else - { - double nVal = n; - for (sal_uLong i = (sal_uLong)k-1; i >= 1; i--) - nVal *= n-(double)i; - PushDouble(nVal); - } - } -} - -void ScInterpreter::ScVariationen2() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScVariationen2" ); - if ( MustHaveParamCount( GetByte(), 2 ) ) - { - double k = ::rtl::math::approxFloor(GetDouble()); - double n = ::rtl::math::approxFloor(GetDouble()); - if (n < 0.0 || k < 0.0 || k > n) - PushIllegalArgument(); - else - PushDouble(pow(n,k)); - } -} - -void ScInterpreter::ScB() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScB" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 3, 4 ) ) - return ; - if (nParamCount == 3) - { - 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 - { - double q = 1.0 - p; - double fFactor = pow(q, n); - if (fFactor == 0.0) - { - fFactor = pow(p, n); - if (fFactor == 0.0) - PushNoValue(); - else - { - sal_uLong max = (sal_uLong) (n - x); - for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) - fFactor *= (n-i)/(i+1)*q/p; - PushDouble(fFactor); - } - } - else - { - sal_uLong max = (sal_uLong) x; - for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) - fFactor *= (n-i)/(i+1)*p/q; - PushDouble(fFactor); - } - } - } - else if (nParamCount == 4) - { - double xe = GetDouble(); - double xs = GetDouble(); - double p = GetDouble(); - double n = GetDouble(); -// alter Stand 300-SC -// if ((xs < n) && (xe < n) && (p < 1.0)) -// { -// double Varianz = sqrt(n * p * (1.0 - p)); -// xs = fabs(xs - (n * p /* / 2.0 STE */ )); -// xe = fabs(xe - (n * p /* / 2.0 STE */ )); -//// STE double nVal = gauss((xs + 0.5) / Varianz) + gauss((xe + 0.5) / Varianz); -// double nVal = fabs(gauss(xs / Varianz) - gauss(xe / Varianz)); -// PushDouble(nVal); -// } - bool bIsValidX = ( 0.0 <= xs && xs <= xe && xe <= n); - if ( bIsValidX && 0.0 < p && p < 1.0 ) - { - double q = 1.0 - p; - double fFactor = pow(q, n); - if (fFactor == 0.0) - { - fFactor = pow(p, n); - if (fFactor == 0.0) - PushNoValue(); - else - { - double fSum = 0.0; - sal_uLong max; - if (xe < (sal_uLong) n) - max = (sal_uLong) (n-xe)-1; - else - max = 0; - sal_uLong i; - for (i = 0; i < max && fFactor > 0.0; i++) - fFactor *= (n-i)/(i+1)*q/p; - if (xs < (sal_uLong) n) - max = (sal_uLong) (n-xs); - else - fSum = fFactor; - for (; i < max && fFactor > 0.0; i++) - { - fFactor *= (n-i)/(i+1)*q/p; - fSum += fFactor; - } - PushDouble(fSum); - } - } - else - { - sal_uLong max; - double fSum; - if ( (sal_uLong) xs == 0) - { - fSum = fFactor; - max = 0; - } - else - { - max = (sal_uLong) xs-1; - fSum = 0.0; - } - sal_uLong i; - for (i = 0; i < max && fFactor > 0.0; i++) - fFactor *= (n-i)/(i+1)*p/q; - if ((sal_uLong)xe == 0) // beide 0 - fSum = fFactor; - else - max = (sal_uLong) xe; - for (; i < max && fFactor > 0.0; i++) - { - fFactor *= (n-i)/(i+1)*p/q; - fSum += fFactor; - } - PushDouble(fSum); - } - } - 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() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScBinomDist" ); - if ( MustHaveParamCount( GetByte(), 4 ) ) - { - double kum = GetDouble(); // 0 oder 1 - double p = GetDouble(); // p - double n = ::rtl::math::approxFloor(GetDouble()); // n - double x = ::rtl::math::approxFloor(GetDouble()); // x - double fFactor, q; - if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0) - PushIllegalArgument(); - else if (kum == 0.0) // Dichte - { - q = 1.0 - p; - fFactor = pow(q, n); - if (fFactor == 0.0) - { - fFactor = pow(p, n); - if (fFactor == 0.0) - PushNoValue(); - else - { - sal_uLong max = (sal_uLong) (n - x); - for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) - fFactor *= (n-i)/(i+1)*q/p; - PushDouble(fFactor); - } - } - else - { - sal_uLong max = (sal_uLong) x; - for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) - fFactor *= (n-i)/(i+1)*p/q; - PushDouble(fFactor); - } - } - else // Verteilung - { - if (n == x) - PushDouble(1.0); - else - { - double fSum; - q = 1.0 - p; - fFactor = pow(q, n); - if (fFactor == 0.0) - { - fFactor = pow(p, n); - if (fFactor == 0.0) - PushNoValue(); - else - { - fSum = 1.0 - fFactor; - sal_uLong max = (sal_uLong) (n - x) - 1; - for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) - { - fFactor *= (n-i)/(i+1)*q/p; - fSum -= fFactor; - } - if (fSum < 0.0) - PushDouble(0.0); - else - PushDouble(fSum); - } - } - else - { - fSum = fFactor; - sal_uLong max = (sal_uLong) x; - for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) - { - fFactor *= (n-i)/(i+1)*p/q; - fSum += fFactor; - } - PushDouble(fSum); - } - } - } - } -} - -void ScInterpreter::ScCritBinom() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScCritBinom" ); - if ( MustHaveParamCount( GetByte(), 3 ) ) - { - double alpha = GetDouble(); // alpha - double p = GetDouble(); // p - double n = ::rtl::math::approxFloor(GetDouble()); - if (n < 0.0 || alpha <= 0.0 || alpha >= 1.0 || p < 0.0 || p > 1.0) - PushIllegalArgument(); - else - { - double q = 1.0 - p; - double fFactor = pow(q,n); - if (fFactor == 0.0) - { - fFactor = pow(p, n); - if (fFactor == 0.0) - PushNoValue(); - else - { - double fSum = 1.0 - fFactor; sal_uLong max = (sal_uLong) n; - sal_uLong i; - - for ( i = 0; i < max && fSum >= alpha; i++) - { - fFactor *= (n-i)/(i+1)*q/p; - fSum -= fFactor; - } - PushDouble(n-i); - } - } - else - { - double fSum = fFactor; sal_uLong max = (sal_uLong) n; - sal_uLong i; - - for ( i = 0; i < max && fSum < alpha; i++) - { - fFactor *= (n-i)/(i+1)*p/q; - fSum += fFactor; - } - PushDouble(i); - } - } - } -} - -void ScInterpreter::ScNegBinomDist() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScNegBinomDist" ); - if ( MustHaveParamCount( GetByte(), 3 ) ) - { - double p = GetDouble(); // p - double r = GetDouble(); // r - double x = GetDouble(); // x - if (r < 0.0 || x < 0.0 || p < 0.0 || p > 1.0) - PushIllegalArgument(); - else - { - double q = 1.0 - p; - double fFactor = pow(p,r); - for (double i = 0.0; i < x; i++) - fFactor *= (i+r)/(i+1.0)*q; - PushDouble(fFactor); - } - } -} - -void ScInterpreter::ScNormDist() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScNormDist" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 3, 4)) - return; - bool bCumulative = nParamCount == 4 ? GetBool() : true; - 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() //expanded, see #i100119# -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScLogNormDist" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 1, 4)) - return; - bool bCumulative = nParamCount == 4 ? GetBool() : true; // 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() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScStdNormDist" ); - PushDouble(integralPhi(GetDouble())); -} - -void ScInterpreter::ScExpDist() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScExpDist" ); - if ( MustHaveParamCount( GetByte(), 3 ) ) - { - double kum = GetDouble(); // 0 oder 1 - double lambda = GetDouble(); // lambda - double x = GetDouble(); // x - if (lambda <= 0.0) - PushIllegalArgument(); - else if (kum == 0.0) // Dichte - { - if (x >= 0.0) - PushDouble(lambda * exp(-lambda*x)); - else - PushInt(0); - } - else // Verteilung - { - if (x > 0.0) - PushDouble(1.0 - exp(-lambda*x)); - else - PushInt(0); - } - } -} - -void ScInterpreter::ScTDist() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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; - } - double R = GetTDist(T, fDF); - if (fFlag == 1.0) - PushDouble(R); - else - PushDouble(2.0*R); -} - -void ScInterpreter::ScFDist() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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::ScChiDist() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScChiDist" ); - 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 ODFF 6.17.10 - { - PushIllegalArgument(); - return; - } - fResult = GetChiDist( fChi, fDF); - if (nGlobalError) - { - PushError( nGlobalError); - return; - } - PushDouble(fResult); -} - -void ScInterpreter::ScWeibull() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScWeibull" ); - if ( MustHaveParamCount( GetByte(), 4 ) ) - { - double kum = GetDouble(); // 0 oder 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) // Dichte - PushDouble(alpha/pow(beta,alpha)*pow(x,alpha-1.0)* - exp(-pow(x/beta,alpha))); - else // Verteilung - PushDouble(1.0 - exp(-pow(x/beta,alpha))); - } -} - -void ScInterpreter::ScPoissonDist() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScPoissonDist" ); - sal_uInt8 nParamCount = GetByte(); - if ( MustHaveParamCount( nParamCount, 2, 3 ) ) - { - bool bCumulative = (nParamCount == 3 ? GetBool() : true); // 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 == 0.0) - PushInt(0); - else - { - if (lambda >712) // 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 == 0.0) - PushInt(1); - else - { - if (lambda > 712 ) // underflow in exp(-lambda) - { // accuracy 12 Digits - PushDouble(GetUpRegIGamma(x+1.0,lambda)); - } - else - { - if (x >= 936.0) // result is always undistinghable from 1 - PushDouble (1.0); - else - { - double fSummand = exp(-lambda); - double fSum = fSummand; - int nEnd = sal::static_int_cast<int>( x ); - for (int i = 1; i <= nEnd; i++) - { - fSummand = (fSummand * lambda)/(double)i; - fSum += fSummand; - } - PushDouble(fSum); - } - } - } - } - } -} - -/** Local function used in the calculation of the hypergeometric distribution. - */ -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. - - 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. - - @author Kohei Yoshida <kohei@openoffice.org> - - @see #i47296# - - */ -void ScInterpreter::ScHypGeomDist() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScHypGeomDist" ); - const size_t nMaxArraySize = 500000; // arbitrary max array size - - if ( !MustHaveParamCount( GetByte(), 4 ) ) - return; - - 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) || (M < x) || (N < n) || (N < M) || (x < n - N + M) ) - { - PushIllegalArgument(); - return; - } - - typedef ::std::vector< double > HypContainer; - HypContainer 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; - } - 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; - } - -#if OSL_DEBUG_LEVEL > 0 - double fCNumLower = N - n - fCNumVarUpper; -#endif - 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() ); - HypContainer::reverse_iterator it1 = cnNumer.rbegin(), it1End = cnNumer.rend(); - HypContainer::reverse_iterator 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; - } - - PushDouble(fFactor); -} - -void ScInterpreter::ScGammaDist() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGammaDist" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 3, 4 ) ) - return; - double bCumulative; - if (nParamCount == 4) - bCumulative = GetBool(); - else - bCumulative = true; - double fBeta = GetDouble(); // scale - double fAlpha = GetDouble(); // shape - double fX = GetDouble(); // x - if (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() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScLogNormInv" ); - if ( MustHaveParamCount( GetByte(), 3 ) ) - { - double sigma = GetDouble(); // Stdabw - double mue = GetDouble(); // Mittelwert - double y = GetDouble(); // y - if (sigma <= 0.0 || y <= 0.0 || y >= 1.0) - PushIllegalArgument(); - else - PushDouble(exp(mue+sigma*gaussinv(y))); - } -} - -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) {} - - double GetValue( double x ) const { return fp - rInt.GetGammaDist(x, fAlpha, fBeta); } -}; - -void ScInterpreter::ScGammaInv() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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(errNoConvergence); - 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) {} - - double GetValue( double x ) const { return fp - rInt.GetBetaDist(x, fAlpha, fBeta); } -}; - -void ScInterpreter::ScBetaInv() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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; - } - if (fP == 0.0) - PushInt(0); - else - { - 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( errNoConvergence); - else - PushDouble(fA + fVal*(fB-fA)); // scale to (A,B) - } -} - - // Achtung: T, F und Chi - // sind monoton fallend, - // deshalb 1-Dist als Funktion - -class ScTDistFunction : public ScDistFunc -{ - ScInterpreter& rInt; - double fp, fDF; - -public: - ScTDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) : - rInt(rI), fp(fpVal), fDF(fDFVal) {} - - double GetValue( double x ) const { return fp - 2 * rInt.GetTDist(x, fDF); } -}; - -void ScInterpreter::ScTInv() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScTInv" ); - if ( !MustHaveParamCount( GetByte(), 2 ) ) - return; - double fDF = ::rtl::math::approxFloor(GetDouble()); - double fP = GetDouble(); - if (fDF < 1.0 || fDF > 1.0E10 || fP <= 0.0 || fP > 1.0 ) - { - PushIllegalArgument(); - return; - } - - bool bConvError; - ScTDistFunction aFunc( *this, fP, fDF ); - double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError ); - if (bConvError) - SetError(errNoConvergence); - PushDouble(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) {} - - double GetValue( double x ) const { return fp - rInt.GetFDist(x, fF1, fF2); } -}; - -void ScInterpreter::ScFInv() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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(errNoConvergence); - 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) {} - - double GetValue( double x ) const { return fp - rInt.GetChiDist(x, fDF); } -}; - -void ScInterpreter::ScChiInv() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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(errNoConvergence); - 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) {} - - double GetValue( double x ) const { 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(errNoConvergence); - PushDouble(fVal); -} - -void ScInterpreter::ScConfidence() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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::ScZTest() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScZTest" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 2, 3 ) ) - return; - double sigma = 0.0, mue, x; - if (nParamCount == 3) - { - sigma = GetDouble(); - if (sigma <= 0.0) - { - PushIllegalArgument(); - return; - } - } - x = GetDouble(); - - double fSum = 0.0; - double fSumSqr = 0.0; - double fVal; - double rValCount = 0.0; - switch (GetStackType()) - { - case formula::svDouble : - { - fVal = GetDouble(); - fSum += fVal; - fSumSqr += fVal*fVal; - rValCount++; - } - break; - case svSingleRef : - { - ScAddress aAdr; - PopSingleRef( aAdr ); - ScBaseCell* pCell = GetCell( aAdr ); - if (HasCellValueData(pCell)) - { - fVal = GetCellValue( aAdr, pCell ); - fSum += fVal; - fSumSqr += fVal*fVal; - rValCount++; - } - } - break; - case svRefList : - case formula::svDoubleRef : - { - short nParam = 1; - size_t nRefInList = 0; - while (nParam-- > 0) - { - ScRange aRange; - sal_uInt16 nErr = 0; - PopDoubleRef( aRange, nParam, nRefInList); - ScValueIterator aValIter(pDok, aRange, glSubTotal); - if (aValIter.GetFirst(fVal, nErr)) - { - fSum += fVal; - fSumSqr += fVal*fVal; - rValCount++; - while ((nErr == 0) && 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->IsString(i)) - { - fVal= pMat->GetDouble(i); - fSum += fVal; - fSumSqr += fVal * fVal; - rValCount++; - } - } - } - } - break; - default : SetError(errIllegalParameter); break; - } - if (rValCount <= 1.0) - PushError( errDivisionByZero); - else - { - mue = fSum/rValCount; - if (nParamCount != 3) - { - sigma = (fSumSqr - fSum*fSum/rValCount)/(rValCount-1.0); - PushDouble(0.5 - gauss((mue-x)/sqrt(sigma/rValCount))); - } - else - PushDouble(0.5 - gauss((mue-x)*sqrt(rValCount)/sigma)); - } -} -bool ScInterpreter::CalculateTest(sal_Bool _bTemplin - ,const SCSIZE nC1, const SCSIZE nC2,const SCSIZE nR1,const SCSIZE nR2 - ,const ScMatrixRef& pMat1,const ScMatrixRef& pMat2 - ,double& fT,double& fF) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateTest" ); - double fCount1 = 0.0; - double fCount2 = 0.0; - double fSum1 = 0.0; - double fSumSqr1 = 0.0; - double fSum2 = 0.0; - double fSumSqr2 = 0.0; - double fVal; - SCSIZE i,j; - for (i = 0; i < nC1; i++) - for (j = 0; j < nR1; j++) - { - if (!pMat1->IsString(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->IsString(i,j)) - { - fVal = pMat2->GetDouble(i,j); - fSum2 += fVal; - fSumSqr2 += fVal * fVal; - fCount2++; - } - } - if (fCount1 < 2.0 || fCount2 < 2.0) - { - PushNoValue(); - return false; - } // if (fCount1 < 2.0 || fCount2 < 2.0) - if ( _bTemplin ) - { - double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0)/fCount1; - double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0)/fCount2; - if (fS1 + fS2 == 0.0) - { - PushNoValue(); - return false; - } - fT = fabs(fSum1/fCount1 - fSum2/fCount2)/sqrt(fS1+fS2); - double c = fS1/(fS1+fS2); - // GetTDist wird mit GetBetaDist berechnet und kommt auch mit nicht ganzzahligen - // Freiheitsgraden klar. Dann stimmt das Ergebnis auch mit Excel ueberein (#52406#): - fF = 1.0/(c*c/(fCount1-1.0)+(1.0-c)*(1.0-c)/(fCount2-1.0)); - } - else - { - // laut Bronstein-Semendjajew - double fS1 = (fSumSqr1 - fSum1*fSum1/fCount1) / (fCount1 - 1.0); // Varianz - double fS2 = (fSumSqr2 - fSum2*fSum2/fCount2) / (fCount2 - 1.0); - fT = fabs( fSum1/fCount1 - fSum2/fCount2 ) / - sqrt( (fCount1-1.0)*fS1 + (fCount2-1.0)*fS2 ) * - sqrt( fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2) ); - fF = fCount1 + fCount2 - 2; - } - return true; -} -void ScInterpreter::ScTTest() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScTTest" ); - if ( !MustHaveParamCount( GetByte(), 4 ) ) - return; - double fTyp = ::rtl::math::approxFloor(GetDouble()); - double fAnz = ::rtl::math::approxFloor(GetDouble()); - if (fAnz != 1.0 && fAnz != 2.0) - { - PushIllegalArgument(); - return; - } - - ScMatrixRef pMat2 = GetMatrix(); - ScMatrixRef pMat1 = GetMatrix(); - if (!pMat1 || !pMat2) - { - PushIllegalParameter(); - return; - } - double fT, fF; - SCSIZE nC1, nC2; - SCSIZE nR1, nR2; - SCSIZE i, j; - pMat1->GetDimensions(nC1, nR1); - pMat2->GetDimensions(nC2, nR2); - if (fTyp == 1.0) - { - if (nC1 != nC2 || nR1 != nR2) - { - PushIllegalArgument(); - return; - } - double fCount = 0.0; - double fSum1 = 0.0; - double fSum2 = 0.0; - double fSumSqrD = 0.0; - double fVal1, fVal2; - for (i = 0; i < nC1; i++) - for (j = 0; j < nR1; j++) - { - if (!pMat1->IsString(i,j) && !pMat2->IsString(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; - } - fT = sqrt(fCount-1.0) * fabs(fSum1 - fSum2) / - sqrt(fCount * fSumSqrD - (fSum1-fSum2)*(fSum1-fSum2)); - fF = fCount - 1.0; - } - else if (fTyp == 2.0) - { - CalculateTest(false,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF); - } - else if (fTyp == 3.0) - { - CalculateTest(sal_True,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF); - } - - else - { - PushIllegalArgument(); - return; - } - if (fAnz == 1.0) - PushDouble(GetTDist(fT, fF)); - else - PushDouble(2.0*GetTDist(fT, fF)); -} - -void ScInterpreter::ScFTest() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScFTest" ); - if ( !MustHaveParamCount( GetByte(), 2 ) ) - return; - ScMatrixRef pMat2 = GetMatrix(); - ScMatrixRef pMat1 = GetMatrix(); - if (!pMat1 || !pMat2) - { - PushIllegalParameter(); - return; - } - SCSIZE nC1, nC2; - SCSIZE nR1, nR2; - SCSIZE i, j; - pMat1->GetDimensions(nC1, nR1); - pMat2->GetDimensions(nC2, nR2); - double fCount1 = 0.0; - double fCount2 = 0.0; - double fSum1 = 0.0; - double fSumSqr1 = 0.0; - double fSum2 = 0.0; - double fSumSqr2 = 0.0; - double fVal; - for (i = 0; i < nC1; i++) - for (j = 0; j < nR1; j++) - { - if (!pMat1->IsString(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->IsString(i,j)) - { - fVal = pMat2->GetDouble(i,j); - fSum2 += fVal; - fSumSqr2 += fVal * fVal; - fCount2++; - } - } - if (fCount1 < 2.0 || fCount2 < 2.0) - { - PushNoValue(); - return; - } - double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0); - double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0); - if (fS1 == 0.0 || fS2 == 0.0) - { - PushNoValue(); - return; - } - double fF, fF1, fF2; - if (fS1 > fS2) - { - fF = fS1/fS2; - fF1 = fCount1-1.0; - fF2 = fCount2-1.0; - } - else - { - fF = fS2/fS1; - fF1 = fCount2-1.0; - fF2 = fCount1-1.0; - } - PushDouble(2.0*GetFDist(fF, fF1, fF2)); -} - -void ScInterpreter::ScChiTest() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScChiTest" ); - if ( !MustHaveParamCount( GetByte(), 2 ) ) - return; - ScMatrixRef pMat2 = GetMatrix(); - ScMatrixRef pMat1 = GetMatrix(); - if (!pMat1 || !pMat2) - { - PushIllegalParameter(); - return; - } - SCSIZE nC1, nC2; - SCSIZE nR1, nR2; - pMat1->GetDimensions(nC1, nR1); - pMat2->GetDimensions(nC2, nR2); - if (nR1 != nR2 || nC1 != nC2) - { - PushIllegalArgument(); - return; - } - double fChi = 0.0; - for (SCSIZE i = 0; i < nC1; i++) - { - for (SCSIZE j = 0; j < nR1; j++) - { - if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j)) - { - double fValX = pMat1->GetDouble(i,j); - double fValE = pMat2->GetDouble(i,j); - fChi += (fValX - fValE) * (fValX - fValE) / fValE; - } - else - { - PushIllegalArgument(); - return; - } - } - } - double fDF; - if (nC1 == 1 || nR1 == 1) - { - fDF = (double)(nC1*nR1 - 1); - if (fDF == 0.0) - { - PushNoValue(); - return; - } - } - else - fDF = (double)(nC1-1)*(double)(nR1-1); - PushDouble(GetChiDist(fChi, fDF)); -} - -void ScInterpreter::ScKurt() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScKurt" ); - double fSum,fCount,vSum; - std::vector<double> values; - if ( !CalculateSkew(fSum,fCount,vSum,values) ) - return; - - if (fCount == 0.0) - { - PushError( errDivisionByZero); - return; - } - - double fMean = fSum / fCount; - - for (size_t i = 0; i < values.size(); i++) - vSum += (values[i] - fMean) * (values[i] - fMean); - - double fStdDev = sqrt(vSum / (fCount - 1.0)); - double dx = 0.0; - double xpower4 = 0.0; - - if (fStdDev == 0.0) - { - PushError( errDivisionByZero); - return; - } - - for (size_t i = 0; i < values.size(); i++) - { - dx = (values[i] - fMean) / fStdDev; - xpower4 = xpower4 + (dx * dx * dx * dx); - } - - double k_d = (fCount - 2.0) * (fCount - 3.0); - double k_l = fCount * (fCount + 1.0) / ((fCount - 1.0) * k_d); - double k_t = 3.0 * (fCount - 1.0) * (fCount - 1.0) / k_d; - - PushDouble(xpower4 * k_l - k_t); -} - -void ScInterpreter::ScHarMean() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScHarMean" ); - short nParamCount = GetByte(); - double nVal = 0.0; - double nValCount = 0.0; - ScAddress aAdr; - ScRange aRange; - size_t nRefInList = 0; - while ((nGlobalError == 0) && (nParamCount-- > 0)) - { - switch (GetStackType()) - { - case formula::svDouble : - { - double x = GetDouble(); - if (x > 0.0) - { - nVal += 1.0/x; - nValCount++; - } - else - SetError( errIllegalArgument); - break; - } - case svSingleRef : - { - PopSingleRef( aAdr ); - ScBaseCell* pCell = GetCell( aAdr ); - if (HasCellValueData(pCell)) - { - double x = GetCellValue( aAdr, pCell ); - if (x > 0.0) - { - nVal += 1.0/x; - nValCount++; - } - else - SetError( errIllegalArgument); - } - break; - } - case formula::svDoubleRef : - case svRefList : - { - sal_uInt16 nErr = 0; - PopDoubleRef( aRange, nParamCount, nRefInList); - double nCellVal; - ScValueIterator aValIter(pDok, aRange, glSubTotal); - if (aValIter.GetFirst(nCellVal, nErr)) - { - if (nCellVal > 0.0) - { - nVal += 1.0/nCellVal; - nValCount++; - } - else - SetError( errIllegalArgument); - SetError(nErr); - while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr)) - { - if (nCellVal > 0.0) - { - nVal += 1.0/nCellVal; - nValCount++; - } - else - SetError( errIllegalArgument); - } - 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( errIllegalArgument); - } - } - else - { - for (SCSIZE nElem = 0; nElem < nCount; nElem++) - if (!pMat->IsString(nElem)) - { - double x = pMat->GetDouble(nElem); - if (x > 0.0) - { - nVal += 1.0/x; - nValCount++; - } - else - SetError( errIllegalArgument); - } - } - } - } - break; - default : SetError(errIllegalParameter); break; - } - } - if (nGlobalError == 0) - PushDouble((double)nValCount/nVal); - else - PushError( nGlobalError); -} - -void ScInterpreter::ScGeoMean() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGeoMean" ); - short nParamCount = GetByte(); - double nVal = 0.0; - double nValCount = 0.0; - ScAddress aAdr; - ScRange aRange; - - size_t nRefInList = 0; - while ((nGlobalError == 0) && (nParamCount-- > 0)) - { - switch (GetStackType()) - { - case formula::svDouble : - { - double x = GetDouble(); - if (x > 0.0) - { - nVal += log(x); - nValCount++; - } - else - SetError( errIllegalArgument); - break; - } - case svSingleRef : - { - PopSingleRef( aAdr ); - ScBaseCell* pCell = GetCell( aAdr ); - if (HasCellValueData(pCell)) - { - double x = GetCellValue( aAdr, pCell ); - if (x > 0.0) - { - nVal += log(x); - nValCount++; - } - else - SetError( errIllegalArgument); - } - break; - } - case formula::svDoubleRef : - case svRefList : - { - sal_uInt16 nErr = 0; - PopDoubleRef( aRange, nParamCount, nRefInList); - double nCellVal; - ScValueIterator aValIter(pDok, aRange, glSubTotal); - if (aValIter.GetFirst(nCellVal, nErr)) - { - if (nCellVal > 0.0) - { - nVal += log(nCellVal); - nValCount++; - } - else - SetError( errIllegalArgument); - SetError(nErr); - while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr)) - { - if (nCellVal > 0.0) - { - nVal += log(nCellVal); - nValCount++; - } - else - SetError( errIllegalArgument); - } - 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 - SetError( errIllegalArgument); - } - } - else - { - for (SCSIZE ui = 0; ui < nCount; ui++) - if (!pMat->IsString(ui)) - { - double x = pMat->GetDouble(ui); - if (x > 0.0) - { - nVal += log(x); - nValCount++; - } - else - SetError( errIllegalArgument); - } - } - } - } - break; - default : SetError(errIllegalParameter); break; - } - } - if (nGlobalError == 0) - PushDouble(exp(nVal / nValCount)); - else - PushError( nGlobalError); -} - -void ScInterpreter::ScStandard() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScStandard" ); - if ( MustHaveParamCount( GetByte(), 3 ) ) - { - double sigma = GetDouble(); - double mue = GetDouble(); - double x = GetDouble(); - if (sigma < 0.0) - PushError( errIllegalArgument); - else if (sigma == 0.0) - PushError( errDivisionByZero); - else - PushDouble((x-mue)/sigma); - } -} -bool ScInterpreter::CalculateSkew(double& fSum,double& fCount,double& vSum,std::vector<double>& values) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateSkew" ); - short nParamCount = GetByte(); - if ( !MustHaveParamCountMin( nParamCount, 1 ) ) - return false; - - fSum = 0.0; - fCount = 0.0; - vSum = 0.0; - double fVal = 0.0; - ScAddress aAdr; - ScRange aRange; - size_t nRefInList = 0; - while (nParamCount-- > 0) - { - switch (GetStackType()) - { - case formula::svDouble : - { - fVal = GetDouble(); - fSum += fVal; - values.push_back(fVal); - fCount++; - } - break; - case svSingleRef : - { - PopSingleRef( aAdr ); - ScBaseCell* pCell = GetCell( aAdr ); - if (HasCellValueData(pCell)) - { - fVal = GetCellValue( aAdr, pCell ); - fSum += fVal; - values.push_back(fVal); - fCount++; - } - } - break; - case formula::svDoubleRef : - case svRefList : - { - PopDoubleRef( aRange, nParamCount, nRefInList); - sal_uInt16 nErr = 0; - ScValueIterator aValIter(pDok, aRange); - if (aValIter.GetFirst(fVal, nErr)) - { - fSum += fVal; - values.push_back(fVal); - fCount++; - SetError(nErr); - while ((nErr == 0) && 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->IsString(nElem)) - { - fVal = pMat->GetDouble(nElem); - fSum += fVal; - values.push_back(fVal); - fCount++; - } - } - } - } - break; - default : - SetError(errIllegalParameter); - break; - } - } - - if (nGlobalError) - { - PushError( nGlobalError); - return false; - } // if (nGlobalError) - return true; -} - -void ScInterpreter::ScSkew() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSkew" ); - double fSum,fCount,vSum; - std::vector<double> values; - if ( !CalculateSkew(fSum,fCount,vSum,values) ) - return; - - double fMean = fSum / fCount; - - for (size_t i = 0; i < values.size(); i++) - vSum += (values[i] - fMean) * (values[i] - fMean); - - double fStdDev = sqrt(vSum / (fCount - 1.0)); - double dx = 0.0; - double xcube = 0.0; - - if (fStdDev == 0) - { - PushIllegalArgument(); - return; - } - - for (size_t i = 0; i < values.size(); i++) - { - dx = (values[i] - fMean) / fStdDev; - xcube = xcube + (dx * dx * dx); - } - - PushDouble(((xcube * fCount) / (fCount - 1.0)) / (fCount - 2.0)); -} - -double ScInterpreter::GetMedian( vector<double> & rArray ) -{ - size_t nSize = rArray.size(); - if (rArray.empty() || nSize == 0 || nGlobalError) - { - SetError( errNoValue); - 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 = rArray.begin() + nMid - 1; - ::std::nth_element( rArray.begin(), iMid, rArray.end()); - return (fUp + *iMid) / 2; - } -} - -void ScInterpreter::ScMedian() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMedian" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCountMin( nParamCount, 1 ) ) - return; - vector<double> aArray; - GetNumberSequenceArray( nParamCount, aArray); - PushDouble( GetMedian( aArray)); -} - -double ScInterpreter::GetPercentile( vector<double> & rArray, double fPercentile ) -{ - size_t nSize = rArray.size(); - if (rArray.empty() || nSize == 0 || nGlobalError) - { - SetError( errNoValue); - return 0.0; - } - - if (nSize == 1) - return rArray[0]; - else - { - size_t nIndex = (size_t)::rtl::math::approxFloor( fPercentile * (nSize-1)); - double fDiff = fPercentile * (nSize-1) - ::rtl::math::approxFloor( fPercentile * (nSize-1)); - OSL_ENSURE(nIndex < nSize, "GetPercentile: wrong index(1)"); - vector<double>::iterator iter = rArray.begin() + nIndex; - ::std::nth_element( rArray.begin(), iter, rArray.end()); - if (fDiff == 0.0) - return *iter; - else - { - OSL_ENSURE(nIndex < nSize-1, "GetPercentile: wrong index(2)"); - double fVal = *iter; - iter = rArray.begin() + nIndex+1; - ::std::nth_element( rArray.begin(), iter, rArray.end()); - return fVal + fDiff * (*iter - fVal); - } - } -} - -void ScInterpreter::ScPercentile() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScPercentile" ); - if ( !MustHaveParamCount( GetByte(), 2 ) ) - return; - double alpha = GetDouble(); - if (alpha < 0.0 || alpha > 1.0) - { - PushIllegalArgument(); - return; - } - vector<double> aArray; - GetNumberSequenceArray( 1, aArray); - PushDouble( GetPercentile( aArray, alpha)); -} - -void ScInterpreter::ScQuartile() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScQuartile" ); - if ( !MustHaveParamCount( GetByte(), 2 ) ) - return; - double fFlag = ::rtl::math::approxFloor(GetDouble()); - if (fFlag < 0.0 || fFlag > 4.0) - { - PushIllegalArgument(); - return; - } - vector<double> aArray; - GetNumberSequenceArray( 1, aArray); - PushDouble( fFlag == 2.0 ? GetMedian( aArray) : GetPercentile( aArray, 0.25 * fFlag)); -} - -void ScInterpreter::ScModalValue() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScModalValue" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCountMin( nParamCount, 1 ) ) - return; - vector<double> aSortArray; - GetSortArray(nParamCount, aSortArray); - SCSIZE nSize = aSortArray.size(); - if (aSortArray.empty() || nSize == 0 || nGlobalError) - 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::CalculateSmallLarge(sal_Bool bSmall) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateSmallLarge" ); - if ( !MustHaveParamCount( GetByte(), 2 ) ) - return; - double f = ::rtl::math::approxFloor(GetDouble()); - if (f < 1.0) - { - PushIllegalArgument(); - return; - } - SCSIZE k = static_cast<SCSIZE>(f); - vector<double> aSortArray; - /* TODO: using nth_element() is best for one single value, but LARGE/SMALL - * actually are defined to return an array of values if an array of - * positions was passed, in which case, depending on the number of values, - * we may or will need a real sorted array again, see #i32345. */ - GetNumberSequenceArray(1, aSortArray); - SCSIZE nSize = aSortArray.size(); - if (aSortArray.empty() || nSize == 0 || nGlobalError || nSize < k) - PushNoValue(); - else - { - // TODO: the sorted case for array: PushDouble( aSortArray[ bSmall ? k-1 : nSize-k ] ); - vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k); - ::std::nth_element( aSortArray.begin(), iPos, aSortArray.end()); - PushDouble( *iPos); - } -} - -void ScInterpreter::ScLarge() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScLarge" ); - CalculateSmallLarge(false); -} - -void ScInterpreter::ScSmall() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSmall" ); - CalculateSmallLarge(sal_True); -} - -void ScInterpreter::ScPercentrank() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScPercentrank" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 2 ) ) - return; - - double fNum = GetDouble(); - vector<double> aSortArray; - GetSortArray(1, aSortArray); - SCSIZE nSize = aSortArray.size(); - if (aSortArray.empty() || nSize == 0 || nGlobalError) - PushNoValue(); - else - { - if (fNum < aSortArray[0] || fNum > aSortArray[nSize-1]) - PushNoValue(); - else if ( nSize == 1 ) - PushDouble(1.0); // fNum == pSortArray[0], see test above - else - { - double fRes; - SCSIZE nOldCount = 0; - double fOldVal = aSortArray[0]; - SCSIZE i; - for (i = 1; i < nSize && aSortArray[i] < fNum; i++) - { - if (aSortArray[i] != fOldVal) - { - nOldCount = i; - fOldVal = aSortArray[i]; - } - } - if (aSortArray[i] != fOldVal) - nOldCount = i; - if (fNum == aSortArray[i]) - fRes = (double)nOldCount/(double)(nSize-1); - else - { - // nOldCount is the count of smaller entries - // fNum is between pSortArray[nOldCount-1] and pSortArray[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 = ( fNum - aSortArray[nOldCount-1] ) / - ( aSortArray[nOldCount] - aSortArray[nOldCount-1] ); - fRes = ( (double)(nOldCount-1)+fFract )/(double)(nSize-1); - } - } - PushDouble(fRes); - } - } -} - -void ScInterpreter::ScTrimMean() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "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); - SCSIZE nSize = aSortArray.size(); - if (aSortArray.empty() || nSize == 0 || nGlobalError) - PushNoValue(); - else - { - sal_uLong nIndex = (sal_uLong) ::rtl::math::approxFloor(alpha*(double)nSize); - if (nIndex % 2 != 0) - nIndex--; - nIndex /= 2; - OSL_ENSURE(nIndex < nSize, "ScTrimMean: falscher Index"); - double fSum = 0.0; - for (SCSIZE i = nIndex; i < nSize-nIndex; i++) - fSum += aSortArray[i]; - PushDouble(fSum/(double)(nSize-2*nIndex)); - } -} - -void ScInterpreter::GetNumberSequenceArray( sal_uInt8 nParamCount, vector<double>& rArray ) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetSortArray" ); - ScAddress aAdr; - ScRange aRange; - short nParam = nParamCount; - size_t nRefInList = 0; - while (nParam-- > 0) - { - switch (GetStackType()) - { - case formula::svDouble : - rArray.push_back( PopDouble()); - break; - case svSingleRef : - { - PopSingleRef( aAdr ); - ScBaseCell* pCell = GetCell( aAdr ); - if (HasCellValueData(pCell)) - rArray.push_back( GetCellValue( aAdr, pCell)); - } - break; - case formula::svDoubleRef : - case svRefList : - { - PopDoubleRef( aRange, nParam, nRefInList); - if (nGlobalError) - break; - - aRange.Justify(); - SCSIZE nCellCount = aRange.aEnd.Col() - aRange.aStart.Col() + 1; - nCellCount *= aRange.aEnd.Row() - aRange.aStart.Row() + 1; - rArray.reserve( rArray.size() + nCellCount); - - sal_uInt16 nErr = 0; - double fCellVal; - ScValueIterator aValIter(pDok, aRange); - if (aValIter.GetFirst( fCellVal, nErr)) - { - rArray.push_back( fCellVal); - SetError(nErr); - while ((nErr == 0) && aValIter.GetNext( fCellVal, nErr)) - rArray.push_back( fCellVal); - SetError(nErr); - } - } - break; - case svMatrix : - case svExternalSingleRef: - case svExternalDoubleRef: - { - ScMatrixRef pMat = GetMatrix(); - if (!pMat) - break; - - SCSIZE nCount = pMat->GetElementCount(); - rArray.reserve( rArray.size() + nCount); - if (pMat->IsNumeric()) - { - for (SCSIZE i = 0; i < nCount; ++i) - rArray.push_back( pMat->GetDouble(i)); - } - else - { - for (SCSIZE i = 0; i < nCount; ++i) - if (!pMat->IsString(i)) - rArray.push_back( pMat->GetDouble(i)); - } - } - break; - default : - PopError(); - SetError( errIllegalParameter); - break; - } - if (nGlobalError) - break; // while - } - // nParam > 0 in case of error, clean stack environment and obtain earlier - // error if there was one. - while (nParam-- > 0) - PopError(); -} - -void ScInterpreter::GetSortArray( sal_uInt8 nParamCount, vector<double>& rSortArray, vector<long>* pIndexOrder ) -{ - GetNumberSequenceArray( nParamCount, rSortArray); - - if (rSortArray.size() > MAX_ANZ_DOUBLE_FOR_SORT) - SetError( errStackOverflow); - else if (rSortArray.empty()) - SetError( errNoValue); - - if (nGlobalError == 0) - QuickSort( rSortArray, pIndexOrder); -} - -static void lcl_QuickSort( long nLo, long nHi, vector<double>& rSortArray, vector<long>* pIndexOrder ) -{ - // If pIndexOrder is not NULL, we assume rSortArray.size() == pIndexOrder->size(). - - using ::std::swap; - - if (nHi - nLo == 1) - { - if (rSortArray[nLo] > rSortArray[nHi]) - { - swap(rSortArray[nLo], rSortArray[nHi]); - if (pIndexOrder) - swap(pIndexOrder->at(nLo), pIndexOrder->at(nHi)); - } - return; - } - - long ni = nLo; - long nj = nHi; - do - { - double fLo = rSortArray[nLo]; - while (ni <= nHi && rSortArray[ni] < fLo) ni++; - while (nj >= nLo && fLo < rSortArray[nj]) nj--; - if (ni <= nj) - { - if (ni != nj) - { - swap(rSortArray[ni], rSortArray[nj]); - if (pIndexOrder) - swap(pIndexOrder->at(ni), pIndexOrder->at(nj)); - } - - ++ni; - --nj; - } - } - while (ni < nj); - - if ((nj - nLo) < (nHi - ni)) - { - if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder); - if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder); - } - else - { - if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder); - if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder); - } -} - -void ScInterpreter::QuickSort( vector<double>& rSortArray, vector<long>* pIndexOrder ) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::QuickSort" ); - long n = static_cast<long>(rSortArray.size()); - - if (pIndexOrder) - { - pIndexOrder->clear(); - pIndexOrder->reserve(n); - for (long i = 0; i < n; ++i) - pIndexOrder->push_back(i); - } - - if (n < 2) - return; - - size_t nValCount = rSortArray.size(); - for (size_t i = 0; (i + 4) <= nValCount-1; i += 4) - { - size_t nInd = rand() % (int) (nValCount-1); - ::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() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRank" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 2, 3 ) ) - return; - sal_Bool bDescending; - if (nParamCount == 3) - bDescending = GetBool(); - else - bDescending = false; - double fCount = 1.0; - sal_Bool bValid = false; - switch (GetStackType()) - { - case formula::svDouble : - { - double x = GetDouble(); - double fVal = GetDouble(); - if (x == fVal) - bValid = sal_True; - break; - } - case svSingleRef : - { - ScAddress aAdr; - PopSingleRef( aAdr ); - double fVal = GetDouble(); - ScBaseCell* pCell = GetCell( aAdr ); - if (HasCellValueData(pCell)) - { - double x = GetCellValue( aAdr, pCell ); - if (x == fVal) - bValid = sal_True; - } - break; - } - case formula::svDoubleRef : - case svRefList : - { - ScRange aRange; - short nParam = 1; - size_t nRefInList = 0; - while (nParam-- > 0) - { - sal_uInt16 nErr = 0; - // Preserve stack until all RefList elements are done! - sal_uInt16 nSaveSP = sp; - PopDoubleRef( aRange, nParam, nRefInList); - if (nParam) - --sp; // simulate pop - double fVal = GetDouble(); - if (nParam) - sp = nSaveSP; - double nCellVal; - ScValueIterator aValIter(pDok, aRange, glSubTotal); - if (aValIter.GetFirst(nCellVal, nErr)) - { - if (nCellVal == fVal) - bValid = sal_True; - else if ((!bDescending && nCellVal > fVal) || - (bDescending && nCellVal < fVal) ) - fCount++; - SetError(nErr); - while ((nErr == 0) && aValIter.GetNext(nCellVal, nErr)) - { - if (nCellVal == fVal) - bValid = sal_True; - else if ((!bDescending && nCellVal > fVal) || - (bDescending && nCellVal < fVal) ) - fCount++; - } - } - SetError(nErr); - } - } - break; - case svMatrix : - case svExternalSingleRef: - case svExternalDoubleRef: - { - ScMatrixRef pMat = GetMatrix(); - double fVal = GetDouble(); - if (pMat) - { - SCSIZE nCount = pMat->GetElementCount(); - if (pMat->IsNumeric()) - { - for (SCSIZE i = 0; i < nCount; i++) - { - double x = pMat->GetDouble(i); - if (x == fVal) - bValid = sal_True; - else if ((!bDescending && x > fVal) || - (bDescending && x < fVal) ) - fCount++; - } - } - else - { - for (SCSIZE i = 0; i < nCount; i++) - if (!pMat->IsString(i)) - { - double x = pMat->GetDouble(i); - if (x == fVal) - bValid = sal_True; - else if ((!bDescending && x > fVal) || - (bDescending && x < fVal) ) - fCount++; - } - } - } - } - break; - default : SetError(errIllegalParameter); break; - } - if (bValid) - PushDouble(fCount); - else - PushNoValue(); -} - -void ScInterpreter::ScAveDev() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScAveDev" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCountMin( nParamCount, 1 ) ) - return; - sal_uInt16 SaveSP = sp; - double nMiddle = 0.0; - double rVal = 0.0; - double rValCount = 0.0; - ScAddress aAdr; - ScRange aRange; - short nParam = nParamCount; - size_t nRefInList = 0; - while (nParam-- > 0) - { - switch (GetStackType()) - { - case formula::svDouble : - rVal += GetDouble(); - rValCount++; - break; - case svSingleRef : - { - PopSingleRef( aAdr ); - ScBaseCell* pCell = GetCell( aAdr ); - if (HasCellValueData(pCell)) - { - rVal += GetCellValue( aAdr, pCell ); - rValCount++; - } - } - break; - case formula::svDoubleRef : - case svRefList : - { - sal_uInt16 nErr = 0; - double nCellVal; - PopDoubleRef( aRange, nParam, nRefInList); - ScValueIterator aValIter(pDok, aRange); - if (aValIter.GetFirst(nCellVal, nErr)) - { - rVal += nCellVal; - rValCount++; - SetError(nErr); - while ((nErr == 0) && 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->IsString(nElem)) - { - rVal += pMat->GetDouble(nElem); - rValCount++; - } - } - } - } - break; - default : - SetError(errIllegalParameter); - break; - } - } - if (nGlobalError) - { - PushError( nGlobalError); - return; - } - nMiddle = rVal / rValCount; - sp = SaveSP; - rVal = 0.0; - nParam = nParamCount; - nRefInList = 0; - while (nParam-- > 0) - { - switch (GetStackType()) - { - case formula::svDouble : - rVal += fabs(GetDouble() - nMiddle); - break; - case svSingleRef : - { - PopSingleRef( aAdr ); - ScBaseCell* pCell = GetCell( aAdr ); - if (HasCellValueData(pCell)) - rVal += fabs(GetCellValue( aAdr, pCell ) - nMiddle); - } - break; - case formula::svDoubleRef : - case svRefList : - { - sal_uInt16 nErr = 0; - double nCellVal; - PopDoubleRef( aRange, nParam, nRefInList); - ScValueIterator aValIter(pDok, aRange); - if (aValIter.GetFirst(nCellVal, nErr)) - { - rVal += (fabs(nCellVal - nMiddle)); - while (aValIter.GetNext(nCellVal, nErr)) - rVal += fabs(nCellVal - nMiddle); - } - } - break; - case svMatrix : - case svExternalSingleRef: - case svExternalDoubleRef: - { - ScMatrixRef pMat = GetMatrix(); - if (pMat) - { - SCSIZE nCount = pMat->GetElementCount(); - if (pMat->IsNumeric()) - { - for (SCSIZE nElem = 0; nElem < nCount; nElem++) - { - rVal += fabs(pMat->GetDouble(nElem) - nMiddle); - } - } - else - { - for (SCSIZE nElem = 0; nElem < nCount; nElem++) - { - if (!pMat->IsString(nElem)) - rVal += fabs(pMat->GetDouble(nElem) - nMiddle); - } - } - } - } - break; - default : SetError(errIllegalParameter); break; - } - } - PushDouble(rVal / rValCount); -} - -void ScInterpreter::ScDevSq() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScDevSq" ); - double nVal; - double nValCount; - GetStVarParams(nVal, nValCount); - PushDouble(nVal); -} - -void ScInterpreter::ScProbability() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScProbability" ); - sal_uInt8 nParamCount = GetByte(); - if ( !MustHaveParamCount( nParamCount, 3, 4 ) ) - return; - double fUp, fLo; - fUp = GetDouble(); - if (nParamCount == 4) - fLo = GetDouble(); - else - fLo = fUp; - if (fLo > fUp) - { - double fTemp = fLo; - fLo = fUp; - fUp = fTemp; - } - ScMatrixRef pMatP = GetMatrix(); - ScMatrixRef pMatW = GetMatrix(); - if (!pMatP || !pMatW) - PushIllegalParameter(); - else - { - SCSIZE nC1, nC2; - SCSIZE nR1, nR2; - pMatP->GetDimensions(nC1, nR1); - pMatW->GetDimensions(nC2, nR2); - if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 || - nC2 == 0 || nR2 == 0) - PushNA(); - else - { - double fSum = 0.0; - double fRes = 0.0; - sal_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( errIllegalArgument); - } - } - if (bStop || fabs(fSum -1.0) > 1.0E-7) - PushNoValue(); - else - PushDouble(fRes); - } - } -} - -void ScInterpreter::ScCorrel() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScCorrel" ); - // This is identical to ScPearson() - ScPearson(); -} - -void ScInterpreter::ScCovar() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScCovar" ); - CalculatePearsonCovar(false,false); -} - -void ScInterpreter::ScPearson() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScPearson" ); - CalculatePearsonCovar(sal_True,false); -} -void ScInterpreter::CalculatePearsonCovar(sal_Bool _bPearson,sal_Bool _bStexy) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculatePearsonCovar" ); - if ( !MustHaveParamCount( GetByte(), 2 ) ) - return; - ScMatrixRef pMat1 = GetMatrix(); - ScMatrixRef pMat2 = GetMatrix(); - if (!pMat1 || !pMat2) - { - PushIllegalParameter(); - return; - } - SCSIZE nC1, nC2; - SCSIZE nR1, nR2; - pMat1->GetDimensions(nC1, nR1); - pMat2->GetDimensions(nC2, nR2); - if (nR1 != nR2 || nC1 != nC2) - { - PushIllegalArgument(); - return; - } - /* #i78250# - * (sum((X-MeanX)(Y-MeanY)))/N equals (SumXY)/N-MeanX*MeanY mathematically, - * but the latter produces wrong results if the absolute values are high, - * for example above 10^8 - */ - double fCount = 0.0; - double fSumX = 0.0; - double fSumY = 0.0; - double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY) - double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2 - double fSumSqrDeltaY = 0.0; // sum of (ValY-MeanY)^2 - for (SCSIZE i = 0; i < nC1; i++) - { - for (SCSIZE j = 0; j < nR1; j++) - { - if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j)) - { - double fValX = pMat1->GetDouble(i,j); - double fValY = pMat2->GetDouble(i,j); - fSumX += fValX; - fSumY += fValY; - fCount++; - } - } - } - if (fCount < (_bStexy ? 3.0 : 1.0)) // fCount==1 is handled by checking denominator later on - PushNoValue(); - else - { - const double fMeanX = fSumX / fCount; - const double fMeanY = fSumY / fCount; - for (SCSIZE i = 0; i < nC1; i++) - { - for (SCSIZE j = 0; j < nR1; j++) - { - if (!pMat1->IsString(i,j) && !pMat2->IsString(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); - } - } - } - } // for (SCSIZE i = 0; i < nC1; i++) - if ( _bPearson ) - { - if (fSumSqrDeltaX == 0.0 || ( !_bStexy && fSumSqrDeltaY == 0.0) ) - PushError( errDivisionByZero); - else if ( _bStexy ) - PushDouble( sqrt( (fSumSqrDeltaY - fSumDeltaXDeltaY * - fSumDeltaXDeltaY / fSumSqrDeltaX) / (fCount-2))); - else - PushDouble( fSumDeltaXDeltaY / sqrt( fSumSqrDeltaX * fSumSqrDeltaY)); - } // if ( _bPearson ) - else - { - PushDouble( fSumDeltaXDeltaY / fCount); - } - } -} - -void ScInterpreter::ScRSQ() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRSQ" ); - // Same as ScPearson()*ScPearson() - ScPearson(); - if (!nGlobalError) - { - switch (GetStackType()) - { - case formula::svDouble: - { - double fVal = PopDouble(); - PushDouble( fVal * fVal); - } - break; - default: - PopError(); - PushNoValue(); - } - } -} - -void ScInterpreter::ScSTEXY() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSTEXY" ); - CalculatePearsonCovar(true,true); -} -void ScInterpreter::CalculateSlopeIntercept(sal_Bool bSlope) -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateSlopeIntercept" ); - if ( !MustHaveParamCount( GetByte(), 2 ) ) - return; - ScMatrixRef pMat1 = GetMatrix(); - ScMatrixRef pMat2 = GetMatrix(); - if (!pMat1 || !pMat2) - { - PushIllegalParameter(); - return; - } - SCSIZE nC1, nC2; - SCSIZE nR1, nR2; - pMat1->GetDimensions(nC1, nR1); - pMat2->GetDimensions(nC2, nR2); - if (nR1 != nR2 || nC1 != nC2) - { - PushIllegalArgument(); - return; - } - // #i78250# numerical stability improved - double fCount = 0.0; - double fSumX = 0.0; - double fSumY = 0.0; - double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY) - double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2 - for (SCSIZE i = 0; i < nC1; i++) - { - for (SCSIZE j = 0; j < nR1; j++) - { - if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j)) - { - double fValX = pMat1->GetDouble(i,j); - double fValY = pMat2->GetDouble(i,j); - fSumX += fValX; - fSumY += fValY; - fCount++; - } - } - } - if (fCount < 1.0) - PushNoValue(); - else - { - double fMeanX = fSumX / fCount; - double fMeanY = fSumY / fCount; - for (SCSIZE i = 0; i < nC1; i++) - { - for (SCSIZE j = 0; j < nR1; j++) - { - if (!pMat1->IsString(i,j) && !pMat2->IsString(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( errDivisionByZero); - else - { - if ( bSlope ) - PushDouble( fSumDeltaXDeltaY / fSumSqrDeltaX); - else - PushDouble( fMeanY - fSumDeltaXDeltaY / fSumSqrDeltaX * fMeanX); - } - } -} - -void ScInterpreter::ScSlope() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSlope" ); - CalculateSlopeIntercept(sal_True); -} - -void ScInterpreter::ScIntercept() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScIntercept" ); - CalculateSlopeIntercept(false); -} - -void ScInterpreter::ScForecast() -{ - RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScForecast" ); - if ( !MustHaveParamCount( GetByte(), 3 ) ) - return; - ScMatrixRef pMat1 = GetMatrix(); - ScMatrixRef pMat2 = GetMatrix(); - if (!pMat1 || !pMat2) - { - PushIllegalParameter(); - return; - } - SCSIZE nC1, nC2; - SCSIZE nR1, nR2; - pMat1->GetDimensions(nC1, nR1); - pMat2->GetDimensions(nC2, nR2); - if (nR1 != nR2 || nC1 != nC2) - { - PushIllegalArgument(); - return; - } - double fVal = GetDouble(); - // #i78250# numerical stability improved - double fCount = 0.0; - double fSumX = 0.0; - double fSumY = 0.0; - double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY) - double fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2 - for (SCSIZE i = 0; i < nC1; i++) - { - for (SCSIZE j = 0; j < nR1; j++) - { - if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j)) - { - double fValX = pMat1->GetDouble(i,j); - double fValY = pMat2->GetDouble(i,j); - fSumX += fValX; - fSumY += fValY; - fCount++; - } - } - } - if (fCount < 1.0) - PushNoValue(); - else - { - double fMeanX = fSumX / fCount; - double fMeanY = fSumY / fCount; - for (SCSIZE i = 0; i < nC1; i++) - { - for (SCSIZE j = 0; j < nR1; j++) - { - if (!pMat1->IsString(i,j) && !pMat2->IsString(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( errDivisionByZero); - else - PushDouble( fMeanY + fSumDeltaXDeltaY / fSumSqrDeltaX * (fVal - fMeanX)); - } -} - -/* vim:set shiftwidth=4 softtabstop=4 expandtab: */ |