horok/libreoffice
1
0
Fork 0
Code Activity
libreoffice/sc/source/core/tool/interpr3.cxx
Daniel Baumann 8e63e14cf6
Adding upstream version 4:25.2.3. 
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
2025-06-22 16:20:04 +02:00

6036 lines
182 KiB
C++
Raw Permalink Blame History

/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
* This file is part of the LibreOffice project.
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* This file incorporates work covered by the following license notice:
*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed
* with this work for additional information regarding copyright
* ownership. The ASF licenses this file to you under the Apache
* License, Version 2.0 (the "License"); you may not use this file
* except in compliance with the License. You may obtain a copy of
* the License at http://www.apache.org/licenses/LICENSE-2.0 .
*/
#include <tools/solar.h>
#include <stdlib.h>
#include <interpre.hxx>
#include <global.hxx>
#include <document.hxx>
#include <dociter.hxx>
#include <matrixoperators.hxx>
#include <scmatrix.hxx>
#include <columniterator.hxx>
#include <unotools/collatorwrapper.hxx>
#include <cassert>
#include <cmath>
#include <memory>
#include <set>
#include <vector>
#include <algorithm>
#include <comphelper/random.hxx>
#include <o3tl/float_int_conversion.hxx>
#include <osl/diagnose.h>
using ::std::vector;
using namespace formula;
/// Two columns of data should be sortable with GetSortArray() and QuickSort()
// This is an arbitrary limit.
static size_t MAX_COUNT_DOUBLE_FOR_SORT(const ScSheetLimits& rSheetLimits)
{
return rSheetLimits.GetMaxRowCount() * 2;
}
const double ScInterpreter::fMaxGammaArgument = 171.624376956302; // found experimental
const double fMachEps = ::std::numeric_limits<double>::epsilon();
namespace {
class ScDistFunc
{
public:
virtual double GetValue(double x) const = 0;
protected:
~ScDistFunc() {}
};
}
// iteration for inverse distributions
//template< class T > double lcl_IterateInverse( const T& rFunction, double x0, double x1, bool& rConvError )
/** u*w<0.0 fails for values near zero */
static bool lcl_HasChangeOfSign( double u, double w )
{
return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);
}
static double lcl_IterateInverse( const ScDistFunc& rFunction, double fAx, double fBx, bool& rConvError )
{
rConvError = false;
const double fYEps = 1.0E-307;
const double fXEps = ::std::numeric_limits<double>::epsilon();
OSL_ENSURE(fAx<fBx, "IterateInverse: wrong interval");
// find enclosing interval
KahanSum fkAx = fAx;
KahanSum fkBx = fBx;
double fAy = rFunction.GetValue(fAx);
double fBy = rFunction.GetValue(fBx);
KahanSum fTemp;
unsigned short nCount;
for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy); nCount++)
{
if (std::abs(fAy) <= std::abs(fBy))
{
fTemp = fkAx;
fkAx += (fkAx - fkBx) * 2.0;
if (fkAx < 0.0)
fkAx = 0.0;
fkBx = fTemp;
fBy = fAy;
fAy = rFunction.GetValue(fkAx.get());
}
else
{
fTemp = fkBx;
fkBx += (fkBx - fkAx) * 2.0;
fkAx = fTemp;
fAy = fBy;
fBy = rFunction.GetValue(fkBx.get());
}
}
fAx = fkAx.get();
fBx = fkBx.get();
if (fAy == 0.0)
return fAx;
if (fBy == 0.0)
return fBx;
if (!lcl_HasChangeOfSign( fAy, fBy))
{
rConvError = true;
return 0.0;
}
// inverse quadric interpolation with additional brackets
// set three points
double fPx = fAx;
double fPy = fAy;
double fQx = fBx;
double fQy = fBy;
double fRx = fAx;
double fRy = fAy;
double fSx = 0.5 * (fAx + fBx); // potential next point
bool bHasToInterpolate = true;
nCount = 0;
while ( nCount < 500 && std::abs(fRy) > fYEps &&
(fBx-fAx) > ::std::max( std::abs(fAx), std::abs(fBx)) * fXEps )
{
if (bHasToInterpolate)
{
if (fPy!=fQy && fQy!=fRy && fRy!=fPy)
{
fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)
+ fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)
+ fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);
bHasToInterpolate = (fAx < fSx) && (fSx < fBx); // inside the brackets?
}
else
bHasToInterpolate = false;
}
if(!bHasToInterpolate)
{
fSx = 0.5 * (fAx + fBx);
// reset points
fQx = fBx; fQy = fBy;
bHasToInterpolate = true;
}
// shift points for next interpolation
fPx = fQx; fQx = fRx; fRx = fSx;
fPy = fQy; fQy = fRy; fRy = rFunction.GetValue(fSx);
// update brackets
if (lcl_HasChangeOfSign( fAy, fRy))
{
fBx = fRx; fBy = fRy;
}
else
{
fAx = fRx; fAy = fRy;
}
// if last iteration brought too small advance, then do bisection next
// time, for safety
bHasToInterpolate = bHasToInterpolate && (std::abs(fRy) * 2.0 <= std::abs(fQy));
++nCount;
}
return fRx;
}
// General functions
void ScInterpreter::ScNoName()
{
PushError(FormulaError::NoName);
}
void ScInterpreter::ScBadName()
{
short nParamCount = GetByte();
while (nParamCount-- > 0)
{
PopError();
}
PushError( FormulaError::NoName);
}
double ScInterpreter::phi(double x)
{
return 0.39894228040143268 * exp(-(x * x) / 2.0);
}
double ScInterpreter::integralPhi(double x)
{ // Using gauss(x)+0.5 has severe cancellation errors for x<-4
return 0.5 * std::erfc(-x * M_SQRT1_2);
}
double ScInterpreter::taylor(const double* pPolynom, sal_uInt16 nMax, double x)
{
KahanSum nVal = pPolynom[nMax];
for (short i = nMax-1; i >= 0; i--)
{
nVal = (nVal * x) + pPolynom[i];
}
return nVal.get();
}
double ScInterpreter::gauss(double x)
{
double xAbs = std::abs(x);
sal_uInt16 xShort = static_cast<sal_uInt16>(::rtl::math::approxFloor(xAbs));
double nVal = 0.0;
if (xShort == 0)
{
static const double t0[] =
{ 0.39894228040143268, -0.06649038006690545, 0.00997355701003582,
-0.00118732821548045, 0.00011543468761616, -0.00000944465625950,
0.00000066596935163, -0.00000004122667415, 0.00000000227352982,
0.00000000011301172, 0.00000000000511243, -0.00000000000021218 };
nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs;
}
else if (xShort <= 2)
{
static const double t2[] =
{ 0.47724986805182079, 0.05399096651318805, -0.05399096651318805,
0.02699548325659403, -0.00449924720943234, -0.00224962360471617,
0.00134977416282970, -0.00011783742691370, -0.00011515930357476,
0.00003704737285544, 0.00000282690796889, -0.00000354513195524,
0.00000037669563126, 0.00000019202407921, -0.00000005226908590,
-0.00000000491799345, 0.00000000366377919, -0.00000000015981997,
-0.00000000017381238, 0.00000000002624031, 0.00000000000560919,
-0.00000000000172127, -0.00000000000008634, 0.00000000000007894 };
nVal = taylor(t2, 23, (xAbs - 2.0));
}
else if (xShort <= 4)
{
static const double t4[] =
{ 0.49996832875816688, 0.00013383022576489, -0.00026766045152977,
0.00033457556441221, -0.00028996548915725, 0.00018178605666397,
-0.00008252863922168, 0.00002551802519049, -0.00000391665839292,
-0.00000074018205222, 0.00000064422023359, -0.00000017370155340,
0.00000000909595465, 0.00000000944943118, -0.00000000329957075,
0.00000000029492075, 0.00000000011874477, -0.00000000004420396,
0.00000000000361422, 0.00000000000143638, -0.00000000000045848 };
nVal = taylor(t4, 20, (xAbs - 4.0));
}
else
{
static const double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };
nVal = 0.5 + phi(xAbs) * taylor(asympt, 4, 1.0 / (xAbs * xAbs)) / xAbs;
}
if (x < 0.0)
return -nVal;
else
return nVal;
}
// #i26836# new gaussinv implementation by Martin Eitzenberger <m.eitzenberger@unix.net>
double ScInterpreter::gaussinv(double x)
{
double q,t,z;
q=x-0.5;
if(std::abs(q)<=.425)
{
t=0.180625-q*q;
z=
q*
(
(
(
(
(
(
(
t*2509.0809287301226727+33430.575583588128105
)
*t+67265.770927008700853
)
*t+45921.953931549871457
)
*t+13731.693765509461125
)
*t+1971.5909503065514427
)
*t+133.14166789178437745
)
*t+3.387132872796366608
)
/
(
(
(
(
(
(
(
t*5226.495278852854561+28729.085735721942674
)
*t+39307.89580009271061
)
*t+21213.794301586595867
)
*t+5394.1960214247511077
)
*t+687.1870074920579083
)
*t+42.313330701600911252
)
*t+1.0
);
}
else
{
if(q>0) t=1-x;
else t=x;
t=sqrt(-log(t));
if(t<=5.0)
{
t+=-1.6;
z=
(
(
(
(
(
(
(
t*7.7454501427834140764e-4+0.0227238449892691845833
)
*t+0.24178072517745061177
)
*t+1.27045825245236838258
)
*t+3.64784832476320460504
)
*t+5.7694972214606914055
)
*t+4.6303378461565452959
)
*t+1.42343711074968357734
)
/
(
(
(
(
(
(
(
t*1.05075007164441684324e-9+5.475938084995344946e-4
)
*t+0.0151986665636164571966
)
*t+0.14810397642748007459
)
*t+0.68976733498510000455
)
*t+1.6763848301838038494
)
*t+2.05319162663775882187
)
*t+1.0
);
}
else
{
t+=-5.0;
z=
(
(
(
(
(
(
(
t*2.01033439929228813265e-7+2.71155556874348757815e-5
)
*t+0.0012426609473880784386
)
*t+0.026532189526576123093
)
*t+0.29656057182850489123
)
*t+1.7848265399172913358
)
*t+5.4637849111641143699
)
*t+6.6579046435011037772
)
/
(
(
(
(
(
(
(
t*2.04426310338993978564e-15+1.4215117583164458887e-7
)
*t+1.8463183175100546818e-5
)
*t+7.868691311456132591e-4
)
*t+0.0148753612908506148525
)
*t+0.13692988092273580531
)
*t+0.59983220655588793769
)
*t+1.0
);
}
if(q<0.0) z=-z;
}
return z;
}
double ScInterpreter::Fakultaet(double x)
{
x = ::rtl::math::approxFloor(x);
if (x < 0.0)
return 0.0;
else if (x == 0.0)
return 1.0;
else if (x <= 170.0)
{
double fTemp = x;
while (fTemp > 2.0)
{
fTemp--;
x *= fTemp;
}
}
else
SetError(FormulaError::NoValue);
return x;
}
double ScInterpreter::BinomKoeff(double n, double k)
{
// this method has been duplicated as BinomialCoefficient()
// in scaddins/source/analysis/analysishelper.cxx
double nVal = 0.0;
k = ::rtl::math::approxFloor(k);
if (n < k)
nVal = 0.0;
else if (k == 0.0)
nVal = 1.0;
else
{
nVal = n/k;
n--;
k--;
while (k > 0.0)
{
nVal *= n/k;
k--;
n--;
}
}
return nVal;
}
// The algorithm is based on lanczos13m53 in lanczos.hpp
// in math library from http://www.boost.org
/** you must ensure fZ>0
Uses a variant of the Lanczos sum with a rational function. */
static double lcl_getLanczosSum(double fZ)
{
static const double fNum[13] ={
23531376880.41075968857200767445163675473,
42919803642.64909876895789904700198885093,
35711959237.35566804944018545154716670596,
17921034426.03720969991975575445893111267,
6039542586.35202800506429164430729792107,
1439720407.311721673663223072794912393972,
248874557.8620541565114603864132294232163,
31426415.58540019438061423162831820536287,
2876370.628935372441225409051620849613599,
186056.2653952234950402949897160456992822,
8071.672002365816210638002902272250613822,
210.8242777515793458725097339207133627117,
2.506628274631000270164908177133837338626
};
static const double fDenom[13] = {
0,
39916800,
120543840,
150917976,
105258076,
45995730,
13339535,
2637558,
357423,
32670,
1925,
66,
1
};
// Horner scheme
double fSumNum;
double fSumDenom;
int nI;
if (fZ<=1.0)
{
fSumNum = fNum[12];
fSumDenom = fDenom[12];
for (nI = 11; nI >= 0; --nI)
{
fSumNum *= fZ;
fSumNum += fNum[nI];
fSumDenom *= fZ;
fSumDenom += fDenom[nI];
}
}
else
// Cancel down with fZ^12; Horner scheme with reverse coefficients
{
double fZInv = 1/fZ;
fSumNum = fNum[0];
fSumDenom = fDenom[0];
for (nI = 1; nI <=12; ++nI)
{
fSumNum *= fZInv;
fSumNum += fNum[nI];
fSumDenom *= fZInv;
fSumDenom += fDenom[nI];
}
}
return fSumNum/fSumDenom;
}
// The algorithm is based on tgamma in gamma.hpp
// in math library from http://www.boost.org
/** You must ensure fZ>0; fZ>171.624376956302 will overflow. */
static double lcl_GetGammaHelper(double fZ)
{
double fGamma = lcl_getLanczosSum(fZ);
const double fg = 6.024680040776729583740234375;
double fZgHelp = fZ + fg - 0.5;
// avoid intermediate overflow
double fHalfpower = pow( fZgHelp, fZ / 2 - 0.25);
fGamma *= fHalfpower;
fGamma /= exp(fZgHelp);
fGamma *= fHalfpower;
if (fZ <= 20.0 && fZ == ::rtl::math::approxFloor(fZ))
fGamma = ::rtl::math::round(fGamma);
return fGamma;
}
// The algorithm is based on tgamma in gamma.hpp
// in math library from http://www.boost.org
/** You must ensure fZ>0 */
static double lcl_GetLogGammaHelper(double fZ)
{
const double fg = 6.024680040776729583740234375;
double fZgHelp = fZ + fg - 0.5;
return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;
}
/** You must ensure non integer arguments for fZ<1 */
double ScInterpreter::GetGamma(double fZ)
{
const double fLogPi = log(M_PI);
const double fLogDblMax = log( ::std::numeric_limits<double>::max());
if (fZ > fMaxGammaArgument)
{
SetError(FormulaError::IllegalFPOperation);
return HUGE_VAL;
}
if (fZ >= 1.0)
return lcl_GetGammaHelper(fZ);
if (fZ >= 0.5) // shift to x>=1 using Gamma(x)=Gamma(x+1)/x
return lcl_GetGammaHelper(fZ+1) / fZ;
if (fZ >= -0.5) // shift to x>=1, might overflow
{
double fLogTest = lcl_GetLogGammaHelper(fZ+2) - std::log1p(fZ) - log( std::abs(fZ));
if (fLogTest >= fLogDblMax)
{
SetError( FormulaError::IllegalFPOperation);
return HUGE_VAL;
}
return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ;
}
// fZ<-0.5
// Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)*sin(pi*x) )
double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( std::abs( ::rtl::math::sin( M_PI*fZ)));
if (fLogDivisor - fLogPi >= fLogDblMax) // underflow
return 0.0;
if (fLogDivisor<0.0)
if (fLogPi - fLogDivisor > fLogDblMax) // overflow
{
SetError(FormulaError::IllegalFPOperation);
return HUGE_VAL;
}
return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( M_PI*fZ) < 0.0) ? -1.0 : 1.0);
}
/** You must ensure fZ>0 */
double ScInterpreter::GetLogGamma(double fZ)
{
if (fZ >= fMaxGammaArgument)
return lcl_GetLogGammaHelper(fZ);
if (fZ >= 1.0)
return log(lcl_GetGammaHelper(fZ));
if (fZ >= 0.5)
return log( lcl_GetGammaHelper(fZ+1) / fZ);
return lcl_GetLogGammaHelper(fZ+2) - std::log1p(fZ) - log(fZ);
}
double ScInterpreter::GetFDist(double x, double fF1, double fF2)
{
double arg = fF2/(fF2+fF1*x);
double alpha = fF2/2.0;
double beta = fF1/2.0;
return GetBetaDist(arg, alpha, beta);
}
double ScInterpreter::GetTDist( double T, double fDF, int nType )
{
switch ( nType )
{
case 1 : // 1-tailed T-distribution
return 0.5 * GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5 );
case 2 : // 2-tailed T-distribution
return GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5);
case 3 : // left-tailed T-distribution (probability density function)
return pow( 1 + ( T * T / fDF ), -( fDF + 1 ) / 2 ) / ( sqrt( fDF ) * GetBeta( 0.5, fDF / 2.0 ) );
case 4 : // left-tailed T-distribution (cumulative distribution function)
double X = fDF / ( T * T + fDF );
double R = 0.5 * GetBetaDist( X, 0.5 * fDF, 0.5 );
return ( T < 0 ? R : 1 - R );
}
SetError( FormulaError::IllegalArgument );
return HUGE_VAL;
}
// for LEGACY.CHIDIST, returns right tail, fDF=degrees of freedom
/** You must ensure fDF>0.0 */
double ScInterpreter::GetChiDist(double fX, double fDF)
{
if (fX <= 0.0)
return 1.0; // see ODFF
else
return GetUpRegIGamma( fDF/2.0, fX/2.0);
}
// ready for ODF 1.2
// for ODF CHISQDIST; cumulative distribution function, fDF=degrees of freedom
// returns left tail
/** You must ensure fDF>0.0 */
double ScInterpreter::GetChiSqDistCDF(double fX, double fDF)
{
if (fX <= 0.0)
return 0.0; // see ODFF
else
return GetLowRegIGamma( fDF/2.0, fX/2.0);
}
double ScInterpreter::GetChiSqDistPDF(double fX, double fDF)
{
// you must ensure fDF is positive integer
double fValue;
if (fX <= 0.0)
return 0.0; // see ODFF
if (fDF*fX > 1391000.0)
{
// intermediate invalid values, use log
fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0) - GetLogGamma(0.5*fDF));
}
else // fDF is small in most cases, we can iterate
{
double fCount;
if (fmod(fDF,2.0)<0.5)
{
// even
fValue = 0.5;
fCount = 2.0;
}
else
{
fValue = 1/sqrt(fX*2*M_PI);
fCount = 1.0;
}
while ( fCount < fDF)
{
fValue *= (fX / fCount);
fCount += 2.0;
}
if (fX>=1425.0) // underflow in e^(-x/2)
fValue = exp(log(fValue)-fX/2);
else
fValue *= exp(-fX/2);
}
return fValue;
}
void ScInterpreter::ScChiSqDist()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
return;
bool bCumulative;
if (nParamCount == 3)
bCumulative = GetBool();
else
bCumulative = true;
double fDF = ::rtl::math::approxFloor(GetDouble());
if (fDF < 1.0)
PushIllegalArgument();
else
{
double fX = GetDouble();
if (bCumulative)
PushDouble(GetChiSqDistCDF(fX,fDF));
else
PushDouble(GetChiSqDistPDF(fX,fDF));
}
}
void ScInterpreter::ScChiSqDist_MS()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 3, 3 ) )
return;
bool bCumulative = GetBool();
double fDF = ::rtl::math::approxFloor( GetDouble() );
if ( fDF < 1.0 || fDF > 1E10 )
PushIllegalArgument();
else
{
double fX = GetDouble();
if ( fX < 0 )
PushIllegalArgument();
else
{
if ( bCumulative )
PushDouble( GetChiSqDistCDF( fX, fDF ) );
else
PushDouble( GetChiSqDistPDF( fX, fDF ) );
}
}
}
void ScInterpreter::ScGamma()
{
double x = GetDouble();
if (x <= 0.0 && x == ::rtl::math::approxFloor(x))
PushIllegalArgument();
else
{
double fResult = GetGamma(x);
if (nGlobalError != FormulaError::NONE)
{
PushError( nGlobalError);
return;
}
PushDouble(fResult);
}
}
void ScInterpreter::ScLogGamma()
{
double x = GetDouble();
if (x > 0.0) // constraint from ODFF
PushDouble( GetLogGamma(x));
else
PushIllegalArgument();
}
double ScInterpreter::GetBeta(double fAlpha, double fBeta)
{
double fA;
double fB;
if (fAlpha > fBeta)
{
fA = fAlpha; fB = fBeta;
}
else
{
fA = fBeta; fB = fAlpha;
}
if (fA+fB < fMaxGammaArgument) // simple case
return GetGamma(fA)/GetGamma(fA+fB)*GetGamma(fB);
// need logarithm
// GetLogGamma is not accurate enough, back to Lanczos for all three
// GetGamma and arrange factors newly.
const double fg = 6.024680040776729583740234375; //see GetGamma
double fgm = fg - 0.5;
double fLanczos = lcl_getLanczosSum(fA);
fLanczos /= lcl_getLanczosSum(fA+fB);
fLanczos *= lcl_getLanczosSum(fB);
double fABgm = fA+fB+fgm;
fLanczos *= sqrt((fABgm/(fA+fgm))/(fB+fgm));
double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
double fTempB = fA/(fB+fgm);
double fResult = exp(-fA * std::log1p(fTempA)
-fB * std::log1p(fTempB)-fgm);
fResult *= fLanczos;
return fResult;
}
// Same as GetBeta but with logarithm
double ScInterpreter::GetLogBeta(double fAlpha, double fBeta)
{
double fA;
double fB;
if (fAlpha > fBeta)
{
fA = fAlpha; fB = fBeta;
}
else
{
fA = fBeta; fB = fAlpha;
}
const double fg = 6.024680040776729583740234375; //see GetGamma
double fgm = fg - 0.5;
double fLanczos = lcl_getLanczosSum(fA);
fLanczos /= lcl_getLanczosSum(fA+fB);
fLanczos *= lcl_getLanczosSum(fB);
double fLogLanczos = log(fLanczos);
double fABgm = fA+fB+fgm;
fLogLanczos += 0.5*(log(fABgm)-log(fA+fgm)-log(fB+fgm));
double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
double fTempB = fA/(fB+fgm);
double fResult = -fA * std::log1p(fTempA)
-fB * std::log1p(fTempB)-fgm;
fResult += fLogLanczos;
return fResult;
}
// beta distribution probability density function
double ScInterpreter::GetBetaDistPDF(double fX, double fA, double fB)
{
// special cases
if (fA == 1.0) // result b*(1-x)^(b-1)
{
if (fB == 1.0)
return 1.0;
if (fB == 2.0)
return -2.0*fX + 2.0;
if (fX == 1.0 && fB < 1.0)
{
SetError(FormulaError::IllegalArgument);
return HUGE_VAL;
}
if (fX <= 0.01)
return fB + fB * std::expm1((fB-1.0) * std::log1p(-fX));
else
return fB * pow(0.5-fX+0.5,fB-1.0);
}
if (fB == 1.0) // result a*x^(a-1)
{
if (fA == 2.0)
return fA * fX;
if (fX == 0.0 && fA < 1.0)
{
SetError(FormulaError::IllegalArgument);
return HUGE_VAL;
}
return fA * pow(fX,fA-1);
}
if (fX <= 0.0)
{
if (fA < 1.0 && fX == 0.0)
{
SetError(FormulaError::IllegalArgument);
return HUGE_VAL;
}
else
return 0.0;
}
if (fX >= 1.0)
{
if (fB < 1.0 && fX == 1.0)
{
SetError(FormulaError::IllegalArgument);
return HUGE_VAL;
}
else
return 0.0;
}
// normal cases; result x^(a-1)*(1-x)^(b-1)/Beta(a,b)
const double fLogDblMax = log( ::std::numeric_limits<double>::max());
const double fLogDblMin = log( ::std::numeric_limits<double>::min());
double fLogY = (fX < 0.1) ? std::log1p(-fX) : log(0.5-fX+0.5);
double fLogX = log(fX);
double fAm1LogX = (fA-1.0) * fLogX;
double fBm1LogY = (fB-1.0) * fLogY;
double fLogBeta = GetLogBeta(fA,fB);
// check whether parts over- or underflow
if ( fAm1LogX < fLogDblMax && fAm1LogX > fLogDblMin
&& fBm1LogY < fLogDblMax && fBm1LogY > fLogDblMin
&& fLogBeta < fLogDblMax && fLogBeta > fLogDblMin
&& fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > fLogDblMin)
return pow(fX,fA-1.0) * pow(0.5-fX+0.5,fB-1.0) / GetBeta(fA,fB);
else // need logarithm;
// might overflow as a whole, but seldom, not worth to pre-detect it
return exp( fAm1LogX + fBm1LogY - fLogBeta);
}
/*
x^a * (1-x)^b
I_x(a,b) = ---------------- * result of ContFrac
a * Beta(a,b)
*/
static double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)
{ // like old version
double a1, b1, a2, b2, fnorm, cfnew, cf;
a1 = 1.0; b1 = 1.0;
b2 = 1.0 - (fA+fB)/(fA+1.0)*fX;
if (b2 == 0.0)
{
a2 = 0.0;
fnorm = 1.0;
cf = 1.0;
}
else
{
a2 = 1.0;
fnorm = 1.0/b2;
cf = a2*fnorm;
}
cfnew = 1.0;
double rm = 1.0;
const double fMaxIter = 50000.0;
// loop security, normal cases converge in less than 100 iterations.
// FIXME: You will get so much iterations for fX near mean,
// I do not know a better algorithm.
bool bfinished = false;
do
{
const double apl2m = fA + 2.0*rm;
const double d2m = rm*(fB-rm)*fX/((apl2m-1.0)*apl2m);
const double d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+1.0));
a1 = (a2+d2m*a1)*fnorm;
b1 = (b2+d2m*b1)*fnorm;
a2 = a1 + d2m1*a2*fnorm;
b2 = b1 + d2m1*b2*fnorm;
if (b2 != 0.0)
{
fnorm = 1.0/b2;
cfnew = a2*fnorm;
bfinished = (std::abs(cf-cfnew) < std::abs(cf)*fMachEps);
}
cf = cfnew;
rm += 1.0;
}
while (rm < fMaxIter && !bfinished);
return cf;
}
// cumulative distribution function, normalized
double ScInterpreter::GetBetaDist(double fXin, double fAlpha, double fBeta)
{
// special cases
if (fXin <= 0.0) // values are valid, see spec
return 0.0;
if (fXin >= 1.0) // values are valid, see spec
return 1.0;
if (fBeta == 1.0)
return pow(fXin, fAlpha);
if (fAlpha == 1.0)
// 1.0 - pow(1.0-fX,fBeta) is not accurate enough
return -std::expm1(fBeta * std::log1p(-fXin));
//FIXME: need special algorithm for fX near fP for large fA,fB
double fResult;
// I use always continued fraction, power series are neither
// faster nor more accurate.
double fY = (0.5-fXin)+0.5;
double flnY = std::log1p(-fXin);
double fX = fXin;
double flnX = log(fXin);
double fA = fAlpha;
double fB = fBeta;
bool bReflect = fXin > fAlpha/(fAlpha+fBeta);
if (bReflect)
{
fA = fBeta;
fB = fAlpha;
fX = fY;
fY = fXin;
flnX = flnY;
flnY = log(fXin);
}
fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);
fResult = fResult/fA;
double fP = fA/(fA+fB);
double fQ = fB/(fA+fB);
double fTemp;
if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97) //found experimental
fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;
else
fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));
fResult *= fTemp;
if (bReflect)
fResult = 0.5 - fResult + 0.5;
if (fResult > 1.0) // ensure valid range
fResult = 1.0;
if (fResult < 0.0)
fResult = 0.0;
return fResult;
}
void ScInterpreter::ScBetaDist()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 3, 6 ) ) // expanded, see #i91547#
return;
double fLowerBound, fUpperBound;
double alpha, beta, x;
bool bIsCumulative;
if (nParamCount == 6)
bIsCumulative = GetBool();
else
bIsCumulative = true;
if (nParamCount >= 5)
fUpperBound = GetDouble();
else
fUpperBound = 1.0;
if (nParamCount >= 4)
fLowerBound = GetDouble();
else
fLowerBound = 0.0;
beta = GetDouble();
alpha = GetDouble();
x = GetDouble();
double fScale = fUpperBound - fLowerBound;
if (fScale <= 0.0 || alpha <= 0.0 || beta <= 0.0)
{
PushIllegalArgument();
return;
}
if (bIsCumulative) // cumulative distribution function
{
// special cases
if (x < fLowerBound)
{
PushDouble(0.0); return; //see spec
}
if (x > fUpperBound)
{
PushDouble(1.0); return; //see spec
}
// normal cases
x = (x-fLowerBound)/fScale; // convert to standard form
PushDouble(GetBetaDist(x, alpha, beta));
return;
}
else // probability density function
{
if (x < fLowerBound || x > fUpperBound)
{
PushDouble(0.0);
return;
}
x = (x-fLowerBound)/fScale;
PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
return;
}
}
/**
Microsoft version has parameters in different order
Also, upper and lowerbound are optional and have default values
and different constraints apply.
Basically, function is identical with ScInterpreter::ScBetaDist()
*/
void ScInterpreter::ScBetaDist_MS()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 4, 6 ) )
return;
double fLowerBound, fUpperBound;
double alpha, beta, x;
bool bIsCumulative;
if (nParamCount == 6)
fUpperBound = GetDouble();
else
fUpperBound = 1.0;
if (nParamCount >= 5)
fLowerBound = GetDouble();
else
fLowerBound = 0.0;
bIsCumulative = GetBool();
beta = GetDouble();
alpha = GetDouble();
x = GetDouble();
if (alpha <= 0.0 || beta <= 0.0 || x < fLowerBound || x > fUpperBound)
{
PushIllegalArgument();
return;
}
double fScale = fUpperBound - fLowerBound;
if (bIsCumulative) // cumulative distribution function
{
x = (x-fLowerBound)/fScale; // convert to standard form
PushDouble(GetBetaDist(x, alpha, beta));
return;
}
else // probability density function
{
x = (x-fLowerBound)/fScale;
PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
return;
}
}
void ScInterpreter::ScPhi()
{
PushDouble(phi(GetDouble()));
}
void ScInterpreter::ScGauss()
{
PushDouble(gauss(GetDouble()));
}
void ScInterpreter::ScFisher()
{
double fVal = GetDouble();
if (std::abs(fVal) >= 1.0)
PushIllegalArgument();
else
PushDouble(::atanh(fVal));
}
void ScInterpreter::ScFisherInv()
{
PushDouble( tanh( GetDouble()));
}
void ScInterpreter::ScFact()
{
double nVal = GetDouble();
if (nVal < 0.0)
PushIllegalArgument();
else
PushDouble(Fakultaet(nVal));
}
void ScInterpreter::ScCombin()
{
if ( MustHaveParamCount( GetByte(), 2 ) )
{
double k = ::rtl::math::approxFloor(GetDouble());
double n = ::rtl::math::approxFloor(GetDouble());
if (k < 0.0 || n < 0.0 || k > n)
PushIllegalArgument();
else
PushDouble(BinomKoeff(n, k));
}
}
void ScInterpreter::ScCombinA()
{
if ( MustHaveParamCount( GetByte(), 2 ) )
{
double k = ::rtl::math::approxFloor(GetDouble());
double n = ::rtl::math::approxFloor(GetDouble());
if (k < 0.0 || n < 0.0 || k > n)
PushIllegalArgument();
else
PushDouble(BinomKoeff(n + k - 1, k));
}
}
void ScInterpreter::ScPermut()
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double k = ::rtl::math::approxFloor(GetDouble());
double n = ::rtl::math::approxFloor(GetDouble());
if (n < 0.0 || k < 0.0 || k > n)
PushIllegalArgument();
else if (k == 0.0)
PushInt(1); // (n! / (n - 0)!) == 1
else
{
double nVal = n;
for (sal_uLong i = static_cast<sal_uLong>(k)-1; i >= 1; i--)
nVal *= n-static_cast<double>(i);
PushDouble(nVal);
}
}
void ScInterpreter::ScPermutationA()
{
if ( MustHaveParamCount( GetByte(), 2 ) )
{
double k = ::rtl::math::approxFloor(GetDouble());
double n = ::rtl::math::approxFloor(GetDouble());
if (n < 0.0 || k < 0.0)
PushIllegalArgument();
else
PushDouble(pow(n,k));
}
}
double ScInterpreter::GetBinomDistPMF(double x, double n, double p)
// used in ScB and ScBinomDist
// preconditions: 0.0 <= x <= n, 0.0 < p < 1.0; x,n integral although double
{
double q = (0.5 - p) + 0.5;
double fFactor = pow(q, n);
if (fFactor <=::std::numeric_limits<double>::min())
{
fFactor = pow(p, n);
if (fFactor <= ::std::numeric_limits<double>::min())
return GetBetaDistPDF(p, x+1.0, n-x+1.0)/(n+1.0);
else
{
sal_uInt32 max = static_cast<sal_uInt32>(n - x);
for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
fFactor *= (n-i)/(i+1)*q/p;
return fFactor;
}
}
else
{
sal_uInt32 max = static_cast<sal_uInt32>(x);
for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
fFactor *= (n-i)/(i+1)*p/q;
return fFactor;
}
}
static double lcl_GetBinomDistRange(double n, double xs,double xe,
double fFactor /* q^n */, double p, double q)
//preconditions: 0.0 <= xs < xe <= n; xs,xe,n integral although double
{
sal_uInt32 i;
// skip summands index 0 to xs-1, start sum with index xs
sal_uInt32 nXs = static_cast<sal_uInt32>( xs );
for (i = 1; i <= nXs && fFactor > 0.0; i++)
fFactor *= (n-i+1)/i * p/q;
KahanSum fSum = fFactor; // Summand xs
sal_uInt32 nXe = static_cast<sal_uInt32>(xe);
for (i = nXs+1; i <= nXe && fFactor > 0.0; i++)
{
fFactor *= (n-i+1)/i * p/q;
fSum += fFactor;
}
return std::min(fSum.get(), 1.0);
}
void ScInterpreter::ScB()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
return ;
if (nParamCount == 3) // mass function
{
double x = ::rtl::math::approxFloor(GetDouble());
double p = GetDouble();
double n = ::rtl::math::approxFloor(GetDouble());
if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
PushIllegalArgument();
else if (p == 0.0)
PushDouble( (x == 0.0) ? 1.0 : 0.0 );
else if ( p == 1.0)
PushDouble( (x == n) ? 1.0 : 0.0);
else
PushDouble(GetBinomDistPMF(x,n,p));
}
else
{ // nParamCount == 4
double xe = ::rtl::math::approxFloor(GetDouble());
double xs = ::rtl::math::approxFloor(GetDouble());
double p = GetDouble();
double n = ::rtl::math::approxFloor(GetDouble());
double q = (0.5 - p) + 0.5;
bool bIsValidX = ( 0.0 <= xs && xs <= xe && xe <= n);
if ( bIsValidX && 0.0 < p && p < 1.0)
{
if (xs == xe) // mass function
PushDouble(GetBinomDistPMF(xs,n,p));
else
{
double fFactor = pow(q, n);
if (fFactor > ::std::numeric_limits<double>::min())
PushDouble(lcl_GetBinomDistRange(n,xs,xe,fFactor,p,q));
else
{
fFactor = pow(p, n);
if (fFactor > ::std::numeric_limits<double>::min())
{
// sum from j=xs to xe {(n choose j) * p^j * q^(n-j)}
// = sum from i = n-xe to n-xs { (n choose i) * q^i * p^(n-i)}
PushDouble(lcl_GetBinomDistRange(n,n-xe,n-xs,fFactor,q,p));
}
else
PushDouble(GetBetaDist(q,n-xe,xe+1.0)-GetBetaDist(q,n-xs+1,xs) );
}
}
}
else
{
if ( bIsValidX ) // not(0<p<1)
{
if ( p == 0.0 )
PushDouble( (xs == 0.0) ? 1.0 : 0.0 );
else if ( p == 1.0 )
PushDouble( (xe == n) ? 1.0 : 0.0 );
else
PushIllegalArgument();
}
else
PushIllegalArgument();
}
}
}
void ScInterpreter::ScBinomDist()
{
if ( !MustHaveParamCount( GetByte(), 4 ) )
return;
bool bIsCum = GetBool(); // false=mass function; true=cumulative
double p = GetDouble();
double n = ::rtl::math::approxFloor(GetDouble());
double x = ::rtl::math::approxFloor(GetDouble());
double q = (0.5 - p) + 0.5; // get one bit more for p near 1.0
if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
{
PushIllegalArgument();
return;
}
if ( p == 0.0)
{
PushDouble( (x==0.0 || bIsCum) ? 1.0 : 0.0 );
return;
}
if ( p == 1.0)
{
PushDouble( (x==n) ? 1.0 : 0.0);
return;
}
if (!bIsCum)
PushDouble( GetBinomDistPMF(x,n,p));
else
{
if (x == n)
PushDouble(1.0);
else
{
double fFactor = pow(q, n);
if (x == 0.0)
PushDouble(fFactor);
else if (fFactor <= ::std::numeric_limits<double>::min())
{
fFactor = pow(p, n);
if (fFactor <= ::std::numeric_limits<double>::min())
PushDouble(GetBetaDist(q,n-x,x+1.0));
else
{
if (fFactor > fMachEps)
{
double fSum = 1.0 - fFactor;
sal_uInt32 max = static_cast<sal_uInt32> (n - x) - 1;
for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
{
fFactor *= (n-i)/(i+1)*q/p;
fSum -= fFactor;
}
PushDouble( (fSum < 0.0) ? 0.0 : fSum );
}
else
PushDouble(lcl_GetBinomDistRange(n,n-x,n,fFactor,q,p));
}
}
else
PushDouble( lcl_GetBinomDistRange(n,0.0,x,fFactor,p,q)) ;
}
}
}
void ScInterpreter::ScCritBinom()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
double alpha = GetDouble();
double p = GetDouble();
double n = ::rtl::math::approxFloor(GetDouble());
if (n < 0.0 || alpha < 0.0 || alpha > 1.0 || p < 0.0 || p > 1.0)
PushIllegalArgument();
else if ( alpha == 0.0 )
PushDouble( 0.0 );
else if ( alpha == 1.0 )
PushDouble( p == 0 ? 0.0 : n );
else
{
double fFactor;
double q = (0.5 - p) + 0.5; // get one bit more for p near 1.0
if ( q > p ) // work from the side where the cumulative curve is
{
// work from 0 upwards
fFactor = pow(q,n);
if (fFactor > ::std::numeric_limits<double>::min())
{
KahanSum fSum = fFactor;
sal_uInt32 max = static_cast<sal_uInt32> (n), i;
for (i = 0; i < max && fSum < alpha; i++)
{
fFactor *= (n-i)/(i+1)*p/q;
fSum += fFactor;
}
PushDouble(i);
}
else
{
// accumulate BinomDist until accumulated BinomDist reaches alpha
KahanSum fSum = 0.0;
sal_uInt32 max = static_cast<sal_uInt32> (n), i;
for (i = 0; i < max && fSum < alpha; i++)
{
const double x = GetBetaDistPDF( p, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
if ( nGlobalError == FormulaError::NONE )
fSum += x;
else
{
PushNoValue();
return;
}
}
assert(i > 0 && "coverity 2023.12.2");
PushDouble( i - 1 );
}
}
else
{
// work from n backwards
fFactor = pow(p, n);
if (fFactor > ::std::numeric_limits<double>::min())
{
KahanSum fSum = 1.0 - fFactor;
sal_uInt32 max = static_cast<sal_uInt32> (n), i;
for (i = 0; i < max && fSum >= alpha; i++)
{
fFactor *= (n-i)/(i+1)*q/p;
fSum -= fFactor;
}
PushDouble(n-i);
}
else
{
// accumulate BinomDist until accumulated BinomDist reaches alpha
KahanSum fSum = 0.0;
sal_uInt32 max = static_cast<sal_uInt32> (n), i;
alpha = 1 - alpha;
for (i = 0; i < max && fSum < alpha; i++)
{
const double x = GetBetaDistPDF( q, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
if ( nGlobalError == FormulaError::NONE )
fSum += x;
else
{
PushNoValue();
return;
}
}
PushDouble( n - i + 1 );
}
}
}
}
void ScInterpreter::ScNegBinomDist()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
double p = GetDouble(); // probability
double s = ::rtl::math::approxFloor(GetDouble()); // No of successes
double f = ::rtl::math::approxFloor(GetDouble()); // No of failures
if ((f + s) <= 1.0 || p < 0.0 || p > 1.0)
PushIllegalArgument();
else
{
double q = 1.0 - p;
double fFactor = pow(p,s);
for (double i = 0.0; i < f; i++)
fFactor *= (i+s)/(i+1.0)*q;
PushDouble(fFactor);
}
}
void ScInterpreter::ScNegBinomDist_MS()
{
if ( !MustHaveParamCount( GetByte(), 4 ) )
return;
bool bCumulative = GetBool();
double p = GetDouble(); // probability
double s = ::rtl::math::approxFloor(GetDouble()); // No of successes
double f = ::rtl::math::approxFloor(GetDouble()); // No of failures
if ( s < 1.0 || f < 0.0 || p < 0.0 || p > 1.0 )
PushIllegalArgument();
else
{
double q = 1.0 - p;
if ( bCumulative )
PushDouble( 1.0 - GetBetaDist( q, f + 1, s ) );
else
{
double fFactor = pow( p, s );
for ( double i = 0.0; i < f; i++ )
fFactor *= ( i + s ) / ( i + 1.0 ) * q;
PushDouble( fFactor );
}
}
}
void ScInterpreter::ScNormDist( int nMinParamCount )
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
return;
bool bCumulative = nParamCount != 4 || GetBool();
double sigma = GetDouble(); // standard deviation
double mue = GetDouble(); // mean
double x = GetDouble(); // x
if (sigma <= 0.0)
{
PushIllegalArgument();
return;
}
if (bCumulative)
PushDouble(integralPhi((x-mue)/sigma));
else
PushDouble(phi((x-mue)/sigma)/sigma);
}
void ScInterpreter::ScLogNormDist( int nMinParamCount ) //expanded, see #i100119# and fdo72158
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
return;
bool bCumulative = nParamCount != 4 || GetBool(); // cumulative
double sigma = nParamCount >= 3 ? GetDouble() : 1.0; // standard deviation
double mue = nParamCount >= 2 ? GetDouble() : 0.0; // mean
double x = GetDouble(); // x
if (sigma <= 0.0)
{
PushIllegalArgument();
return;
}
if (bCumulative)
{ // cumulative
if (x <= 0.0)
PushDouble(0.0);
else
PushDouble(integralPhi((log(x)-mue)/sigma));
}
else
{ // density
if (x <= 0.0)
PushIllegalArgument();
else
PushDouble(phi((log(x)-mue)/sigma)/sigma/x);
}
}
void ScInterpreter::ScStdNormDist()
{
PushDouble(integralPhi(GetDouble()));
}
void ScInterpreter::ScStdNormDist_MS()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 2 ) )
return;
bool bCumulative = GetBool(); // cumulative
double x = GetDouble(); // x
if ( bCumulative )
PushDouble( integralPhi( x ) );
else
PushDouble( exp( - pow( x, 2 ) / 2 ) / sqrt( 2 * M_PI ) );
}
void ScInterpreter::ScExpDist()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
double kum = GetDouble(); // 0 or 1
double lambda = GetDouble(); // lambda
double x = GetDouble(); // x
if (lambda <= 0.0)
PushIllegalArgument();
else if (kum == 0.0) // density
{
if (x >= 0.0)
PushDouble(lambda * exp(-lambda*x));
else
PushInt(0);
}
else // distribution
{
if (x > 0.0)
PushDouble(1.0 - exp(-lambda*x));
else
PushInt(0);
}
}
void ScInterpreter::ScTDist()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
double fFlag = ::rtl::math::approxFloor(GetDouble());
double fDF = ::rtl::math::approxFloor(GetDouble());
double T = GetDouble();
if (fDF < 1.0 || T < 0.0 || (fFlag != 1.0 && fFlag != 2.0) )
{
PushIllegalArgument();
return;
}
PushDouble( GetTDist( T, fDF, static_cast<int>(fFlag) ) );
}
void ScInterpreter::ScTDist_T( int nTails )
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double fDF = ::rtl::math::approxFloor( GetDouble() );
double fT = GetDouble();
if ( fDF < 1.0 || ( nTails == 2 && fT < 0.0 ) )
{
PushIllegalArgument();
return;
}
double fRes = GetTDist( fT, fDF, nTails );
if ( nTails == 1 && fT < 0.0 )
PushDouble( 1.0 - fRes ); // tdf#105937, right tail, negative X
else
PushDouble( fRes );
}
void ScInterpreter::ScTDist_MS()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
bool bCumulative = GetBool();
double fDF = ::rtl::math::approxFloor( GetDouble() );
double T = GetDouble();
if ( fDF < 1.0 )
{
PushIllegalArgument();
return;
}
PushDouble( GetTDist( T, fDF, ( bCumulative ? 4 : 3 ) ) );
}
void ScInterpreter::ScFDist()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
double fF2 = ::rtl::math::approxFloor(GetDouble());
double fF1 = ::rtl::math::approxFloor(GetDouble());
double fF = GetDouble();
if (fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10)
{
PushIllegalArgument();
return;
}
PushDouble(GetFDist(fF, fF1, fF2));
}
void ScInterpreter::ScFDist_LT()
{
int nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
return;
bool bCum;
if ( nParamCount == 3 )
bCum = true;
else if ( IsMissing() )
{
bCum = true;
Pop();
}
else
bCum = GetBool();
double fF2 = ::rtl::math::approxFloor( GetDouble() );
double fF1 = ::rtl::math::approxFloor( GetDouble() );
double fF = GetDouble();
if ( fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 )
{
PushIllegalArgument();
return;
}
if ( bCum )
{
// left tail cumulative distribution
PushDouble( 1.0 - GetFDist( fF, fF1, fF2 ) );
}
else
{
// probability density function
PushDouble( pow( fF1 / fF2, fF1 / 2 ) * pow( fF, ( fF1 / 2 ) - 1 ) /
( pow( ( 1 + ( fF * fF1 / fF2 ) ), ( fF1 + fF2 ) / 2 ) *
GetBeta( fF1 / 2, fF2 / 2 ) ) );
}
}
void ScInterpreter::ScChiDist( bool bODFF )
{
double fResult;
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double fDF = ::rtl::math::approxFloor(GetDouble());
double fChi = GetDouble();
if ( fDF < 1.0 // x<=0 returns 1, see ODFF1.2 6.18.11
|| ( !bODFF && fChi < 0 ) ) // Excel does not accept negative fChi
{
PushIllegalArgument();
return;
}
fResult = GetChiDist( fChi, fDF);
if (nGlobalError != FormulaError::NONE)
{
PushError( nGlobalError);
return;
}
PushDouble(fResult);
}
void ScInterpreter::ScWeibull()
{
if ( !MustHaveParamCount( GetByte(), 4 ) )
return;
double kum = GetDouble(); // 0 or 1
double beta = GetDouble(); // beta
double alpha = GetDouble(); // alpha
double x = GetDouble(); // x
if (alpha <= 0.0 || beta <= 0.0 || x < 0.0)
PushIllegalArgument();
else if (kum == 0.0) // Density
PushDouble(alpha/pow(beta,alpha)*pow(x,alpha-1.0)*
exp(-pow(x/beta,alpha)));
else // Distribution
PushDouble(1.0 - exp(-pow(x/beta,alpha)));
}
void ScInterpreter::ScPoissonDist( bool bODFF )
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, ( bODFF ? 2 : 3 ), 3 ) )
return;
bool bCumulative = nParamCount != 3 || GetBool(); // default cumulative
double lambda = GetDouble(); // Mean
double x = ::rtl::math::approxFloor(GetDouble()); // discrete distribution
if (lambda <= 0.0 || x < 0.0)
PushIllegalArgument();
else if (!bCumulative) // Probability mass function
{
if (lambda >712.0) // underflow in exp(-lambda)
{ // accuracy 11 Digits
PushDouble( exp(x*log(lambda)-lambda-GetLogGamma(x+1.0)));
}
else
{
double fPoissonVar = 1.0;
for ( double f = 0.0; f < x; ++f )
fPoissonVar *= lambda / ( f + 1.0 );
PushDouble( fPoissonVar * exp( -lambda ) );
}
}
else // Cumulative distribution function
{
if (lambda > 712.0) // underflow in exp(-lambda)
{ // accuracy 12 Digits
PushDouble(GetUpRegIGamma(x+1.0,lambda));
}
else
{
if (x >= 936.0) // result is always indistinguishable from 1
PushDouble (1.0);
else
{
double fSummand = std::exp(-lambda);
KahanSum fSum = fSummand;
int nEnd = sal::static_int_cast<int>( x );
for (int i = 1; i <= nEnd; i++)
{
fSummand = (fSummand * lambda)/static_cast<double>(i);
fSum += fSummand;
}
PushDouble(fSum.get());
}
}
}
}
/** Local function used in the calculation of the hypergeometric distribution.
*/
static void lcl_PutFactorialElements( ::std::vector< double >& cn, double fLower, double fUpper, double fBase )
{
for ( double i = fLower; i <= fUpper; ++i )
{
double fVal = fBase - i;
if ( fVal > 1.0 )
cn.push_back( fVal );
}
}
/** Calculates a value of the hypergeometric distribution.
@see #i47296#
This function has an extra argument bCumulative,
which only calculates the non-cumulative distribution and
which is optional in Calc and mandatory with Excel's HYPGEOM.DIST()
@see fdo#71722
@see tdf#102948, make Calc function ODFF1.2-compliant
@see tdf#117041, implement note at bottom of ODFF1.2 par.6.18.37
*/
void ScInterpreter::ScHypGeomDist( int nMinParamCount )
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, nMinParamCount, 5 ) )
return;
bool bCumulative = ( nParamCount == 5 && GetBool() );
double N = ::rtl::math::approxFloor(GetDouble());
double M = ::rtl::math::approxFloor(GetDouble());
double n = ::rtl::math::approxFloor(GetDouble());
double x = ::rtl::math::approxFloor(GetDouble());
if ( (x < 0.0) || (n < x) || (N < n) || (N < M) || (M < 0.0) )
{
PushIllegalArgument();
return;
}
KahanSum fVal = 0.0;
for ( int i = ( bCumulative ? 0 : x ); i <= x && nGlobalError == FormulaError::NONE; i++ )
{
if ( (n - i <= N - M) && (i <= M) )
fVal += GetHypGeomDist( i, n, M, N );
}
PushDouble( fVal.get() );
}
/** Calculates a value of the hypergeometric distribution.
The algorithm is designed to avoid unnecessary multiplications and division
by expanding all factorial elements (9 of them). It is done by excluding
those ranges that overlap in the numerator and the denominator. This allows
for a fast calculation for large values which would otherwise cause an overflow
in the intermediate values.
@see #i47296#
*/
double ScInterpreter::GetHypGeomDist( double x, double n, double M, double N )
{
const size_t nMaxArraySize = 500000; // arbitrary max array size
std::vector<double> cnNumer, cnDenom;
size_t nEstContainerSize = static_cast<size_t>( x + ::std::min( n, M ) );
size_t nMaxSize = ::std::min( cnNumer.max_size(), nMaxArraySize );
if ( nEstContainerSize > nMaxSize )
{
PushNoValue();
return 0;
}
cnNumer.reserve( nEstContainerSize + 10 );
cnDenom.reserve( nEstContainerSize + 10 );
// Trim coefficient C first
double fCNumVarUpper = N - n - M + x - 1.0;
double fCDenomVarLower = 1.0;
if ( N - n - M + x >= M - x + 1.0 )
{
fCNumVarUpper = M - x - 1.0;
fCDenomVarLower = N - n - 2.0*(M - x) + 1.0;
}
double fCNumLower = N - n - fCNumVarUpper;
double fCDenomUpper = N - n - M + x + 1.0 - fCDenomVarLower;
double fDNumVarLower = n - M;
if ( n >= M + 1.0 )
{
if ( N - M < n + 1.0 )
{
// Case 1
if ( N - n < n + 1.0 )
{
// no overlap
lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
lcl_PutFactorialElements( cnDenom, 0.0, N - n - 1.0, N );
}
else
{
// overlap
OSL_ENSURE( fCNumLower < n + 1.0, "ScHypGeomDist: wrong assertion" );
lcl_PutFactorialElements( cnNumer, N - 2.0*n, fCNumVarUpper, N - n );
lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
}
OSL_ENSURE( fCDenomUpper <= N - M, "ScHypGeomDist: wrong assertion" );
if ( fCDenomUpper < n - x + 1.0 )
// no overlap
lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
else
{
// overlap
lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
fCDenomUpper = n - x;
fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
}
}
else
{
// Case 2
if ( n > M - 1.0 )
{
// no overlap
lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
}
else
{
lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
}
OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
if ( fCDenomUpper < n - x + 1.0 )
// no overlap
lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - n + x, N - M + 1.0 );
else
{
lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
fCDenomUpper = n - x;
fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
}
}
OSL_ENSURE( fCDenomUpper <= M, "ScHypGeomDist: wrong assertion" );
}
else
{
if ( N - M < M + 1.0 )
{
// Case 3
if ( N - n < M + 1.0 )
{
// No overlap
lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
lcl_PutFactorialElements( cnDenom, 0.0, N - M - 1.0, N );
}
else
{
lcl_PutFactorialElements( cnNumer, N - n - M, fCNumVarUpper, N - n );
lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
}
if ( n - x + 1.0 > fCDenomUpper )
// No overlap
lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
else
{
// Overlap
lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );
fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
fCDenomUpper = n - x;
}
}
else
{
// Case 4
OSL_ENSURE( M >= n - x, "ScHypGeomDist: wrong assertion" );
OSL_ENSURE( M - x <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
if ( N - n < N - M + 1.0 )
{
// No overlap
lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
}
else
{
// Overlap
OSL_ENSURE( fCNumLower <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
}
if ( n - x + 1.0 > fCDenomUpper )
// No overlap
lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - n + x, N - M + 1.0 );
else if ( M >= fCDenomUpper )
{
lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
fCDenomUpper = n - x;
fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
}
else
{
OSL_ENSURE( M <= fCDenomUpper, "ScHypGeomDist: wrong assertion" );
lcl_PutFactorialElements( cnDenom, fCDenomVarLower, N - n - 2.0*M + x,
N - n - M + x + 1.0 );
fCDenomUpper = n - x;
fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
}
}
OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );
fDNumVarLower = 0.0;
}
double nDNumVarUpper = fCDenomUpper < x + 1.0 ? n - x - 1.0 : n - fCDenomUpper - 1.0;
double nDDenomVarLower = fCDenomUpper < x + 1.0 ? fCDenomVarLower : N - n - M + 1.0;
lcl_PutFactorialElements( cnNumer, fDNumVarLower, nDNumVarUpper, n );
lcl_PutFactorialElements( cnDenom, nDDenomVarLower, N - n - M + x, N - n - M + x + 1.0 );
::std::sort( cnNumer.begin(), cnNumer.end() );
::std::sort( cnDenom.begin(), cnDenom.end() );
auto it1 = cnNumer.rbegin(), it1End = cnNumer.rend();
auto it2 = cnDenom.rbegin(), it2End = cnDenom.rend();
double fFactor = 1.0;
for ( ; it1 != it1End || it2 != it2End; )
{
double fEnum = 1.0, fDenom = 1.0;
if ( it1 != it1End )
fEnum = *it1++;
if ( it2 != it2End )
fDenom = *it2++;
fFactor *= fEnum / fDenom;
}
return fFactor;
}
void ScInterpreter::ScGammaDist( bool bODFF )
{
sal_uInt8 nMinParamCount = ( bODFF ? 3 : 4 );
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
return;
bool bCumulative;
if (nParamCount == 4)
bCumulative = GetBool();
else
bCumulative = true;
double fBeta = GetDouble(); // scale
double fAlpha = GetDouble(); // shape
double fX = GetDouble(); // x
if ((!bODFF && fX < 0) || fAlpha <= 0.0 || fBeta <= 0.0)
PushIllegalArgument();
else
{
if (bCumulative) // distribution
PushDouble( GetGammaDist( fX, fAlpha, fBeta));
else // density
PushDouble( GetGammaDistPDF( fX, fAlpha, fBeta));
}
}
void ScInterpreter::ScNormInv()
{
if ( MustHaveParamCount( GetByte(), 3 ) )
{
double sigma = GetDouble();
double mue = GetDouble();
double x = GetDouble();
if (sigma <= 0.0 || x < 0.0 || x > 1.0)
PushIllegalArgument();
else if (x == 0.0 || x == 1.0)
PushNoValue();
else
PushDouble(gaussinv(x)*sigma + mue);
}
}
void ScInterpreter::ScSNormInv()
{
double x = GetDouble();
if (x < 0.0 || x > 1.0)
PushIllegalArgument();
else if (x == 0.0 || x == 1.0)
PushNoValue();
else
PushDouble(gaussinv(x));
}
void ScInterpreter::ScLogNormInv()
{
sal_uInt8 nParamCount = GetByte();
if ( MustHaveParamCount( nParamCount, 1, 3 ) )
{
double fSigma = ( nParamCount == 3 ? GetDouble() : 1.0 ); // Stddev
double fMue = ( nParamCount >= 2 ? GetDouble() : 0.0 ); // Mean
double fP = GetDouble(); // p
if ( fSigma <= 0.0 || fP <= 0.0 || fP >= 1.0 )
PushIllegalArgument();
else
PushDouble( exp( fMue + fSigma * gaussinv( fP ) ) );
}
}
class ScGammaDistFunction : public ScDistFunc
{
ScInterpreter& rInt;
double fp, fAlpha, fBeta;
public:
ScGammaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
virtual ~ScGammaDistFunction() {}
double GetValue( double x ) const override { return fp - rInt.GetGammaDist(x, fAlpha, fBeta); }
};
void ScInterpreter::ScGammaInv()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
double fBeta = GetDouble();
double fAlpha = GetDouble();
double fP = GetDouble();
if (fAlpha <= 0.0 || fBeta <= 0.0 || fP < 0.0 || fP >= 1.0 )
{
PushIllegalArgument();
return;
}
if (fP == 0.0)
PushInt(0);
else
{
bool bConvError;
ScGammaDistFunction aFunc( *this, fP, fAlpha, fBeta );
double fStart = fAlpha * fBeta;
double fVal = lcl_IterateInverse( aFunc, fStart*0.5, fStart, bConvError );
if (bConvError)
SetError(FormulaError::NoConvergence);
PushDouble(fVal);
}
}
class ScBetaDistFunction : public ScDistFunc
{
ScInterpreter& rInt;
double fp, fAlpha, fBeta;
public:
ScBetaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}
virtual ~ScBetaDistFunction() {}
double GetValue( double x ) const override { return fp - rInt.GetBetaDist(x, fAlpha, fBeta); }
};
void ScInterpreter::ScBetaInv()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 3, 5 ) )
return;
double fP, fA, fB, fAlpha, fBeta;
if (nParamCount == 5)
fB = GetDouble();
else
fB = 1.0;
if (nParamCount >= 4)
fA = GetDouble();
else
fA = 0.0;
fBeta = GetDouble();
fAlpha = GetDouble();
fP = GetDouble();
if (fP < 0.0 || fP > 1.0 || fA >= fB || fAlpha <= 0.0 || fBeta <= 0.0)
{
PushIllegalArgument();
return;
}
bool bConvError;
ScBetaDistFunction aFunc( *this, fP, fAlpha, fBeta );
// 0..1 as range for iteration so it isn't extended beyond the valid range
double fVal = lcl_IterateInverse( aFunc, 0.0, 1.0, bConvError );
if (bConvError)
PushError( FormulaError::NoConvergence);
else
PushDouble(fA + fVal*(fB-fA)); // scale to (A,B)
}
// Note: T, F, and Chi are
// monotonically decreasing,
// therefore 1-Dist as function
class ScTDistFunction : public ScDistFunc
{
ScInterpreter& rInt;
double fp, fDF;
int nT;
public:
ScTDistFunction( ScInterpreter& rI, double fpVal, double fDFVal, int nType ) :
rInt( rI ), fp( fpVal ), fDF( fDFVal ), nT( nType ) {}
virtual ~ScTDistFunction() {}
double GetValue( double x ) const override { return fp - rInt.GetTDist( x, fDF, nT ); }
};
void ScInterpreter::ScTInv( int nType )
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double fDF = ::rtl::math::approxFloor(GetDouble());
double fP = GetDouble();
if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
{
PushIllegalArgument();
return;
}
if ( nType == 4 ) // left-tailed cumulative t-distribution
{
if ( fP == 1.0 )
PushIllegalArgument();
else if ( fP < 0.5 )
PushDouble( -GetTInv( 1 - fP, fDF, nType ) );
else
PushDouble( GetTInv( fP, fDF, nType ) );
}
else
PushDouble( GetTInv( fP, fDF, nType ) );
};
double ScInterpreter::GetTInv( double fAlpha, double fSize, int nType )
{
bool bConvError;
ScTDistFunction aFunc( *this, fAlpha, fSize, nType );
double fVal = lcl_IterateInverse( aFunc, fSize * 0.5, fSize, bConvError );
if (bConvError)
SetError(FormulaError::NoConvergence);
return fVal;
}
class ScFDistFunction : public ScDistFunc
{
ScInterpreter& rInt;
double fp, fF1, fF2;
public:
ScFDistFunction( ScInterpreter& rI, double fpVal, double fF1Val, double fF2Val ) :
rInt(rI), fp(fpVal), fF1(fF1Val), fF2(fF2Val) {}
virtual ~ScFDistFunction() {}
double GetValue( double x ) const override { return fp - rInt.GetFDist(x, fF1, fF2); }
};
void ScInterpreter::ScFInv()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
double fF2 = ::rtl::math::approxFloor(GetDouble());
double fF1 = ::rtl::math::approxFloor(GetDouble());
double fP = GetDouble();
if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
{
PushIllegalArgument();
return;
}
bool bConvError;
ScFDistFunction aFunc( *this, fP, fF1, fF2 );
double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
if (bConvError)
SetError(FormulaError::NoConvergence);
PushDouble(fVal);
}
void ScInterpreter::ScFInv_LT()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
double fF2 = ::rtl::math::approxFloor(GetDouble());
double fF1 = ::rtl::math::approxFloor(GetDouble());
double fP = GetDouble();
if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
{
PushIllegalArgument();
return;
}
bool bConvError;
ScFDistFunction aFunc( *this, ( 1.0 - fP ), fF1, fF2 );
double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
if (bConvError)
SetError(FormulaError::NoConvergence);
PushDouble(fVal);
}
class ScChiDistFunction : public ScDistFunc
{
ScInterpreter& rInt;
double fp, fDF;
public:
ScChiDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
rInt(rI), fp(fpVal), fDF(fDFVal) {}
virtual ~ScChiDistFunction() {}
double GetValue( double x ) const override { return fp - rInt.GetChiDist(x, fDF); }
};
void ScInterpreter::ScChiInv()
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double fDF = ::rtl::math::approxFloor(GetDouble());
double fP = GetDouble();
if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
{
PushIllegalArgument();
return;
}
bool bConvError;
ScChiDistFunction aFunc( *this, fP, fDF );
double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
if (bConvError)
SetError(FormulaError::NoConvergence);
PushDouble(fVal);
}
/***********************************************/
class ScChiSqDistFunction : public ScDistFunc
{
ScInterpreter& rInt;
double fp, fDF;
public:
ScChiSqDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
rInt(rI), fp(fpVal), fDF(fDFVal) {}
virtual ~ScChiSqDistFunction() {}
double GetValue( double x ) const override { return fp - rInt.GetChiSqDistCDF(x, fDF); }
};
void ScInterpreter::ScChiSqInv()
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double fDF = ::rtl::math::approxFloor(GetDouble());
double fP = GetDouble();
if (fDF < 1.0 || fP < 0.0 || fP >= 1.0 )
{
PushIllegalArgument();
return;
}
bool bConvError;
ScChiSqDistFunction aFunc( *this, fP, fDF );
double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
if (bConvError)
SetError(FormulaError::NoConvergence);
PushDouble(fVal);
}
void ScInterpreter::ScConfidence()
{
if ( MustHaveParamCount( GetByte(), 3 ) )
{
double n = ::rtl::math::approxFloor(GetDouble());
double sigma = GetDouble();
double alpha = GetDouble();
if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
PushIllegalArgument();
else
PushDouble( gaussinv(1.0-alpha/2.0) * sigma/sqrt(n) );
}
}
void ScInterpreter::ScConfidenceT()
{
if ( MustHaveParamCount( GetByte(), 3 ) )
{
double n = ::rtl::math::approxFloor(GetDouble());
double sigma = GetDouble();
double alpha = GetDouble();
if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
PushIllegalArgument();
else if (n == 1.0) // for interoperability with Excel
PushError(FormulaError::DivisionByZero);
else
PushDouble( sigma * GetTInv( alpha, n - 1, 2 ) / sqrt( n ) );
}
}
void ScInterpreter::ScZTest()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
return;
double sigma = 0.0, x;
if (nParamCount == 3)
{
sigma = GetDouble();
if (sigma <= 0.0)
{
PushIllegalArgument();
return;
}
}
x = GetDouble();
KahanSum fSum = 0.0;
KahanSum fSumSqr = 0.0;
double fVal;
double rValCount = 0.0;
switch (GetStackType())
{
case svDouble :
{
fVal = GetDouble();
fSum += fVal;
fSumSqr += fVal*fVal;
rValCount++;
}
break;
case svSingleRef :
{
ScAddress aAdr;
PopSingleRef( aAdr );
ScRefCellValue aCell(mrDoc, aAdr);
if (aCell.hasNumeric())
{
fVal = GetCellValue(aAdr, aCell);
fSum += fVal;
fSumSqr += fVal*fVal;
rValCount++;
}
}
break;
case svRefList :
case svDoubleRef :
{
short nParam = 1;
size_t nRefInList = 0;
while (nParam-- > 0)
{
ScRange aRange;
FormulaError nErr = FormulaError::NONE;
PopDoubleRef( aRange, nParam, nRefInList);
ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
if (aValIter.GetFirst(fVal, nErr))
{
fSum += fVal;
fSumSqr += fVal*fVal;
rValCount++;
while ((nErr == FormulaError::NONE) && aValIter.GetNext(fVal, nErr))
{
fSum += fVal;
fSumSqr += fVal*fVal;
rValCount++;
}
SetError(nErr);
}
}
}
break;
case svMatrix :
case svExternalSingleRef:
case svExternalDoubleRef:
{
ScMatrixRef pMat = GetMatrix();
if (pMat)
{
SCSIZE nCount = pMat->GetElementCount();
if (pMat->IsNumeric())
{
for ( SCSIZE i = 0; i < nCount; i++ )
{
fVal= pMat->GetDouble(i);
fSum += fVal;
fSumSqr += fVal * fVal;
rValCount++;
}
}
else
{
for (SCSIZE i = 0; i < nCount; i++)
if (!pMat->IsStringOrEmpty(i))
{
fVal= pMat->GetDouble(i);
fSum += fVal;
fSumSqr += fVal * fVal;
rValCount++;
}
}
}
}
break;
default : SetError(FormulaError::IllegalParameter); break;
}
if (rValCount <= 1.0)
PushError( FormulaError::DivisionByZero);
else
{
double mue = fSum.get()/rValCount;
if (nParamCount != 3)
{
sigma = (fSumSqr - fSum*fSum/rValCount).get()/(rValCount-1.0);
if (sigma == 0.0)
{
PushError(FormulaError::DivisionByZero);
return;
}
PushDouble(0.5 - gauss((mue-x)/sqrt(sigma/rValCount)));
}
else
PushDouble(0.5 - gauss((mue-x)*sqrt(rValCount)/sigma));
}
}
bool ScInterpreter::CalculateTest(bool _bTemplin
,const SCSIZE nC1, const SCSIZE nC2,const SCSIZE nR1,const SCSIZE nR2
,const ScMatrixRef& pMat1,const ScMatrixRef& pMat2
,double& fT,double& fF)
{
double fCount1 = 0.0;
double fCount2 = 0.0;
KahanSum fSum1 = 0.0;
KahanSum fSumSqr1 = 0.0;
KahanSum fSum2 = 0.0;
KahanSum fSumSqr2 = 0.0;
double fVal;
SCSIZE i,j;
for (i = 0; i < nC1; i++)
for (j = 0; j < nR1; j++)
{
if (!pMat1->IsStringOrEmpty(i,j))
{
fVal = pMat1->GetDouble(i,j);
fSum1 += fVal;
fSumSqr1 += fVal * fVal;
fCount1++;
}
}
for (i = 0; i < nC2; i++)
for (j = 0; j < nR2; j++)
{
if (!pMat2->IsStringOrEmpty(i,j))
{
fVal = pMat2->GetDouble(i,j);
fSum2 += fVal;
fSumSqr2 += fVal * fVal;
fCount2++;
}
}
if (fCount1 < 2.0 || fCount2 < 2.0)
{
PushNoValue();
return false;
} // if (fCount1 < 2.0 || fCount2 < 2.0)
if ( _bTemplin )
{
double fS1 = (fSumSqr1-fSum1*fSum1/fCount1).get() / (fCount1-1.0) / fCount1;
double fS2 = (fSumSqr2-fSum2*fSum2/fCount2).get() / (fCount2-1.0) / fCount2;
if (fS1 + fS2 == 0.0)
{
PushNoValue();
return false;
}
fT = std::abs(( fSum1/fCount1 - fSum2/fCount2 ).get())/sqrt(fS1+fS2);
double c = fS1/(fS1+fS2);
// GetTDist is calculated via GetBetaDist and also works with non-integral
// degrees of freedom. The result matches Excel
fF = 1.0/(c*c/(fCount1-1.0)+(1.0-c)*(1.0-c)/(fCount2-1.0));
}
else
{
// according to Bronstein-Semendjajew
double fS1 = (fSumSqr1 - fSum1*fSum1/fCount1).get() / (fCount1 - 1.0); // Variance
double fS2 = (fSumSqr2 - fSum2*fSum2/fCount2).get() / (fCount2 - 1.0);
fT = std::abs( fSum1.get()/fCount1 - fSum2.get()/fCount2 ) /
sqrt( (fCount1-1.0)*fS1 + (fCount2-1.0)*fS2 ) *
sqrt( fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2) );
fF = fCount1 + fCount2 - 2;
}
return true;
}
void ScInterpreter::ScTTest()
{
if ( !MustHaveParamCount( GetByte(), 4 ) )
return;
double fTyp = ::rtl::math::approxFloor(GetDouble());
double fTails = ::rtl::math::approxFloor(GetDouble());
if (fTails != 1.0 && fTails != 2.0)
{
PushIllegalArgument();
return;
}
ScMatrixRef pMat2 = GetMatrix();
ScMatrixRef pMat1 = GetMatrix();
if (!pMat1 || !pMat2)
{
PushIllegalParameter();
return;
}
double fT, fF;
SCSIZE nC1, nC2;
SCSIZE nR1, nR2;
SCSIZE i, j;
pMat1->GetDimensions(nC1, nR1);
pMat2->GetDimensions(nC2, nR2);
if (fTyp == 1.0)
{
if (nC1 != nC2 || nR1 != nR2)
{
PushIllegalArgument();
return;
}
double fCount = 0.0;
KahanSum fSum1 = 0.0;
KahanSum fSum2 = 0.0;
KahanSum fSumSqrD = 0.0;
double fVal1, fVal2;
for (i = 0; i < nC1; i++)
for (j = 0; j < nR1; j++)
{
if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
{
fVal1 = pMat1->GetDouble(i,j);
fVal2 = pMat2->GetDouble(i,j);
fSum1 += fVal1;
fSum2 += fVal2;
fSumSqrD += (fVal1 - fVal2)*(fVal1 - fVal2);
fCount++;
}
}
if (fCount < 1.0)
{
PushNoValue();
return;
}
KahanSum fSumD = fSum1 - fSum2;
double fDivider = ( fSumSqrD*fCount - fSumD*fSumD ).get();
if ( fDivider == 0.0 )
{
PushError(FormulaError::DivisionByZero);
return;
}
fT = std::abs(fSumD.get()) * sqrt((fCount-1.0) / fDivider);
fF = fCount - 1.0;
}
else if (fTyp == 2.0)
{
if (!CalculateTest(false,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF))
return; // error was pushed
}
else if (fTyp == 3.0)
{
if (!CalculateTest(true,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF))
return; // error was pushed
}
else
{
PushIllegalArgument();
return;
}
PushDouble( GetTDist( fT, fF, static_cast<int>(fTails) ) );
}
void ScInterpreter::ScFTest()
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
ScMatrixRef pMat2 = GetMatrix();
ScMatrixRef pMat1 = GetMatrix();
if (!pMat1 || !pMat2)
{
PushIllegalParameter();
return;
}
auto aVal1 = pMat1->CollectKahan(sc::op::kOpSumAndSumSquare);
auto aVal2 = pMat2->CollectKahan(sc::op::kOpSumAndSumSquare);
double fCount1 = aVal1.mnCount;
double fCount2 = aVal2.mnCount;
KahanSum fSum1 = aVal1.maAccumulator[0];
KahanSum fSumSqr1 = aVal1.maAccumulator[1];
KahanSum fSum2 = aVal2.maAccumulator[0];
KahanSum fSumSqr2 = aVal2.maAccumulator[1];
if (fCount1 < 2.0 || fCount2 < 2.0)
{
PushNoValue();
return;
}
double fS1 = (fSumSqr1-fSum1*fSum1/fCount1).get() / (fCount1-1.0);
double fS2 = (fSumSqr2-fSum2*fSum2/fCount2).get() / (fCount2-1.0);
if (fS1 == 0.0 || fS2 == 0.0)
{
PushNoValue();
return;
}
double fF, fF1, fF2;
if (fS1 > fS2)
{
fF = fS1/fS2;
fF1 = fCount1-1.0;
fF2 = fCount2-1.0;
}
else
{
fF = fS2/fS1;
fF1 = fCount2-1.0;
fF2 = fCount1-1.0;
}
double fFcdf = GetFDist(fF, fF1, fF2);
PushDouble(2.0*std::min(fFcdf, 1.0 - fFcdf));
}
void ScInterpreter::ScChiTest()
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
ScMatrixRef pMat2 = GetMatrix();
ScMatrixRef pMat1 = GetMatrix();
if (!pMat1 || !pMat2)
{
PushIllegalParameter();
return;
}
SCSIZE nC1, nC2;
SCSIZE nR1, nR2;
pMat1->GetDimensions(nC1, nR1);
pMat2->GetDimensions(nC2, nR2);
if (nR1 != nR2 || nC1 != nC2)
{
PushIllegalArgument();
return;
}
KahanSum fChi = 0.0;
bool bEmpty = true;
for (SCSIZE i = 0; i < nC1; i++)
{
for (SCSIZE j = 0; j < nR1; j++)
{
if (!(pMat1->IsEmpty(i,j) || pMat2->IsEmpty(i,j)))
{
bEmpty = false;
if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
{
double fValX = pMat1->GetDouble(i,j);
double fValE = pMat2->GetDouble(i,j);
if ( fValE == 0.0 )
{
PushError(FormulaError::DivisionByZero);
return;
}
// These fTemp values guard against a failure when compiled
// with optimization (using g++ 4.8.2 on tinderbox 71-TDF),
// where ((fValX - fValE) * (fValX - fValE)) with
// fValE==1e+308 should had produced Infinity but did
// not, instead the result of divide() then was 1e+308.
volatile double fTemp1 = (fValX - fValE) * (fValX - fValE);
double fTemp2 = fTemp1;
if (std::isinf(fTemp2))
{
PushError(FormulaError::NoConvergence);
return;
}
fChi += sc::divide( fTemp2, fValE);
}
else
{
PushIllegalArgument();
return;
}
}
}
}
if ( bEmpty )
{
// not in ODFF1.2, but for interoperability with Excel
PushIllegalArgument();
return;
}
double fDF;
if (nC1 == 1 || nR1 == 1)
{
fDF = static_cast<double>(nC1*nR1 - 1);
if (fDF == 0.0)
{
PushNoValue();
return;
}
}
else
fDF = static_cast<double>(nC1-1)*static_cast<double>(nR1-1);
PushDouble(GetChiDist(fChi.get(), fDF));
}
void ScInterpreter::ScKurt()
{
KahanSum fSum;
double fCount;
std::vector<double> values;
if ( !CalculateSkew(fSum, fCount, values) )
return;
// ODF 1.2 constraints: # of numbers >= 4
if (fCount < 4.0)
{
// for interoperability with Excel
PushError( FormulaError::DivisionByZero);
return;
}
KahanSum vSum;
double fMean = fSum.get() / fCount;
for (double v : values)
vSum += (v - fMean) * (v - fMean);
double fStdDev = sqrt(vSum.get() / (fCount - 1.0));
if (fStdDev == 0.0)
{
PushError( FormulaError::DivisionByZero);
return;
}
KahanSum xpower4 = 0.0;
for (double v : values)
{
double dx = (v - fMean) / fStdDev;
xpower4 += (dx * dx) * (dx * dx);
}
double k_d = (fCount - 2.0) * (fCount - 3.0);
double k_l = fCount * (fCount + 1.0) / ((fCount - 1.0) * k_d);
double k_t = 3.0 * (fCount - 1.0) * (fCount - 1.0) / k_d;
PushDouble(xpower4.get() * k_l - k_t);
}
void ScInterpreter::ScHarMean()
{
short nParamCount = GetByte();
KahanSum nVal = 0.0;
double nValCount = 0.0;
ScAddress aAdr;
ScRange aRange;
size_t nRefInList = 0;
while ((nGlobalError == FormulaError::NONE) && (nParamCount-- > 0))
{
switch (GetStackType())
{
case svDouble :
{
double x = GetDouble();
if (x > 0.0)
{
nVal += 1.0/x;
nValCount++;
}
else
SetError( FormulaError::IllegalArgument);
break;
}
case svSingleRef :
{
PopSingleRef( aAdr );
ScRefCellValue aCell(mrDoc, aAdr);
if (aCell.hasNumeric())
{
double x = GetCellValue(aAdr, aCell);
if (x > 0.0)
{
nVal += 1.0/x;
nValCount++;
}
else
SetError( FormulaError::IllegalArgument);
}
break;
}
case svDoubleRef :
case svRefList :
{
FormulaError nErr = FormulaError::NONE;
PopDoubleRef( aRange, nParamCount, nRefInList);
double nCellVal;
ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
if (aValIter.GetFirst(nCellVal, nErr))
{
if (nCellVal > 0.0)
{
nVal += 1.0/nCellVal;
nValCount++;
}
else
SetError( FormulaError::IllegalArgument);
SetError(nErr);
while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
{
if (nCellVal > 0.0)
{
nVal += 1.0/nCellVal;
nValCount++;
}
else
SetError( FormulaError::IllegalArgument);
}
SetError(nErr);
}
}
break;
case svMatrix :
case svExternalSingleRef:
case svExternalDoubleRef:
{
ScMatrixRef pMat = GetMatrix();
if (pMat)
{
SCSIZE nCount = pMat->GetElementCount();
if (pMat->IsNumeric())
{
for (SCSIZE nElem = 0; nElem < nCount; nElem++)
{
double x = pMat->GetDouble(nElem);
if (x > 0.0)
{
nVal += 1.0/x;
nValCount++;
}
else
SetError( FormulaError::IllegalArgument);
}
}
else
{
for (SCSIZE nElem = 0; nElem < nCount; nElem++)
if (!pMat->IsStringOrEmpty(nElem))
{
double x = pMat->GetDouble(nElem);
if (x > 0.0)
{
nVal += 1.0/x;
nValCount++;
}
else
SetError( FormulaError::IllegalArgument);
}
}
}
}
break;
default : SetError(FormulaError::IllegalParameter); break;
}
}
if (nGlobalError == FormulaError::NONE)
PushDouble( nValCount / nVal.get() );
else
PushError( nGlobalError);
}
void ScInterpreter::ScGeoMean()
{
short nParamCount = GetByte();
KahanSum nVal = 0.0;
double nValCount = 0.0;
ScAddress aAdr;
ScRange aRange;
size_t nRefInList = 0;
while ((nGlobalError == FormulaError::NONE) && (nParamCount-- > 0))
{
switch (GetStackType())
{
case svDouble :
{
double x = GetDouble();
if (x > 0.0)
{
nVal += log(x);
nValCount++;
}
else if ( x == 0.0 )
{
// value of 0 means that function result will be 0
while ( nParamCount-- > 0 )
PopError();
PushDouble( 0.0 );
return;
}
else
SetError( FormulaError::IllegalArgument);
break;
}
case svSingleRef :
{
PopSingleRef( aAdr );
ScRefCellValue aCell(mrDoc, aAdr);
if (aCell.hasNumeric())
{
double x = GetCellValue(aAdr, aCell);
if (x > 0.0)
{
nVal += log(x);
nValCount++;
}
else if ( x == 0.0 )
{
// value of 0 means that function result will be 0
while ( nParamCount-- > 0 )
PopError();
PushDouble( 0.0 );
return;
}
else
SetError( FormulaError::IllegalArgument);
}
break;
}
case svDoubleRef :
case svRefList :
{
FormulaError nErr = FormulaError::NONE;
PopDoubleRef( aRange, nParamCount, nRefInList);
double nCellVal;
ScValueIterator aValIter(mrContext, aRange, mnSubTotalFlags);
if (aValIter.GetFirst(nCellVal, nErr))
{
if (nCellVal > 0.0)
{
nVal += log(nCellVal);
nValCount++;
}
else if ( nCellVal == 0.0 )
{
// value of 0 means that function result will be 0
while ( nParamCount-- > 0 )
PopError();
PushDouble( 0.0 );
return;
}
else
SetError( FormulaError::IllegalArgument);
SetError(nErr);
while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
{
if (nCellVal > 0.0)
{
nVal += log(nCellVal);
nValCount++;
}
else if ( nCellVal == 0.0 )
{
// value of 0 means that function result will be 0
while ( nParamCount-- > 0 )
PopError();
PushDouble( 0.0 );
return;
}
else
SetError( FormulaError::IllegalArgument);
}
SetError(nErr);
}
}
break;
case svMatrix :
case svExternalSingleRef:
case svExternalDoubleRef:
{
ScMatrixRef pMat = GetMatrix();
if (pMat)
{
SCSIZE nCount = pMat->GetElementCount();
if (pMat->IsNumeric())
{
for (SCSIZE ui = 0; ui < nCount; ui++)
{
double x = pMat->GetDouble(ui);
if (x > 0.0)
{
nVal += log(x);
nValCount++;
}
else if ( x == 0.0 )
{
// value of 0 means that function result will be 0
while ( nParamCount-- > 0 )
PopError();
PushDouble( 0.0 );
return;
}
else
SetError( FormulaError::IllegalArgument);
}
}
else
{
for (SCSIZE ui = 0; ui < nCount; ui++)
{
if (!pMat->IsStringOrEmpty(ui))
{
double x = pMat->GetDouble(ui);
if (x > 0.0)
{
nVal += log(x);
nValCount++;
}
else if ( x == 0.0 )
{
// value of 0 means that function result will be 0
while ( nParamCount-- > 0 )
PopError();
PushDouble( 0.0 );
return;
}
else
SetError( FormulaError::IllegalArgument);
}
}
}
}
}
break;
default : SetError(FormulaError::IllegalParameter); break;
}
}
if (nGlobalError == FormulaError::NONE)
PushDouble(exp(nVal.get() / nValCount));
else
PushError( nGlobalError);
}
void ScInterpreter::ScStandard()
{
if ( MustHaveParamCount( GetByte(), 3 ) )
{
double sigma = GetDouble();
double mue = GetDouble();
double x = GetDouble();
if (sigma < 0.0)
PushError( FormulaError::IllegalArgument);
else if (sigma == 0.0)
PushError( FormulaError::DivisionByZero);
else
PushDouble((x-mue)/sigma);
}
}
bool ScInterpreter::CalculateSkew(KahanSum& fSum, double& fCount, std::vector<double>& values)
{
short nParamCount = GetByte();
if ( !MustHaveParamCountMin( nParamCount, 1 ) )
return false;
fSum = 0.0;
fCount = 0.0;
double fVal = 0.0;
ScAddress aAdr;
ScRange aRange;
size_t nRefInList = 0;
while (nParamCount-- > 0)
{
switch (GetStackType())
{
case svDouble :
{
fVal = GetDouble();
fSum += fVal;
values.push_back(fVal);
fCount++;
}
break;
case svSingleRef :
{
PopSingleRef( aAdr );
ScRefCellValue aCell(mrDoc, aAdr);
if (aCell.hasNumeric())
{
fVal = GetCellValue(aAdr, aCell);
fSum += fVal;
values.push_back(fVal);
fCount++;
}
}
break;
case svDoubleRef :
case svRefList :
{
PopDoubleRef( aRange, nParamCount, nRefInList);
FormulaError nErr = FormulaError::NONE;
ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
if (aValIter.GetFirst(fVal, nErr))
{
fSum += fVal;
values.push_back(fVal);
fCount++;
SetError(nErr);
while ((nErr == FormulaError::NONE) && aValIter.GetNext(fVal, nErr))
{
fSum += fVal;
values.push_back(fVal);
fCount++;
}
SetError(nErr);
}
}
break;
case svMatrix :
case svExternalSingleRef:
case svExternalDoubleRef:
{
ScMatrixRef pMat = GetMatrix();
if (pMat)
{
SCSIZE nCount = pMat->GetElementCount();
if (pMat->IsNumeric())
{
for (SCSIZE nElem = 0; nElem < nCount; nElem++)
{
fVal = pMat->GetDouble(nElem);
fSum += fVal;
values.push_back(fVal);
fCount++;
}
}
else
{
for (SCSIZE nElem = 0; nElem < nCount; nElem++)
if (!pMat->IsStringOrEmpty(nElem))
{
fVal = pMat->GetDouble(nElem);
fSum += fVal;
values.push_back(fVal);
fCount++;
}
}
}
}
break;
default :
SetError(FormulaError::IllegalParameter);
break;
}
}
if (nGlobalError != FormulaError::NONE)
{
PushError( nGlobalError);
return false;
} // if (nGlobalError != FormulaError::NONE)
return true;
}
void ScInterpreter::CalculateSkewOrSkewp( bool bSkewp )
{
KahanSum fSum;
double fCount;
std::vector<double> values;
if (!CalculateSkew( fSum, fCount, values))
return;
// SKEW/SKEWP's constraints: they require at least three numbers
if (fCount < 3.0)
{
// for interoperability with Excel
PushError(FormulaError::DivisionByZero);
return;
}
KahanSum vSum;
double fMean = fSum.get() / fCount;
for (double v : values)
vSum += (v - fMean) * (v - fMean);
double fStdDev = sqrt( vSum.get() / (bSkewp ? fCount : (fCount - 1.0)));
if (fStdDev == 0)
{
PushIllegalArgument();
return;
}
KahanSum xcube = 0.0;
for (double v : values)
{
double dx = (v - fMean) / fStdDev;
xcube += dx * dx * dx;
}
if (bSkewp)
PushDouble( xcube.get() / fCount );
else
PushDouble( ((xcube.get() * fCount) / (fCount - 1.0)) / (fCount - 2.0) );
}
void ScInterpreter::ScSkew()
{
CalculateSkewOrSkewp( false );
}
void ScInterpreter::ScSkewp()
{
CalculateSkewOrSkewp( true );
}
double ScInterpreter::GetMedian( vector<double> & rArray )
{
size_t nSize = rArray.size();
if (nSize == 0 || nGlobalError != FormulaError::NONE)
{
SetError( FormulaError::NoValue);
return 0.0;
}
// Upper median.
size_t nMid = nSize / 2;
vector<double>::iterator iMid = rArray.begin() + nMid;
::std::nth_element( rArray.begin(), iMid, rArray.end());
if (nSize & 1)
return *iMid; // Lower and upper median are equal.
else
{
double fUp = *iMid;
// Lower median.
iMid = ::std::max_element( rArray.begin(), rArray.begin() + nMid);
return (fUp + *iMid) / 2;
}
}
void ScInterpreter::ScMedian()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCountMin( nParamCount, 1 ) )
return;
vector<double> aArray;
GetNumberSequenceArray( nParamCount, aArray, false );
PushDouble( GetMedian( aArray));
}
double ScInterpreter::GetPercentile( vector<double> & rArray, double fPercentile )
{
size_t nSize = rArray.size();
if (nSize == 1)
return rArray[0];
else
{
size_t nIndex = static_cast<size_t>(::rtl::math::approxFloor( fPercentile * (nSize-1)));
double fDiff = fPercentile * (nSize-1) - ::rtl::math::approxFloor( fPercentile * (nSize-1));
OSL_ENSURE(nIndex < nSize, "GetPercentile: wrong index(1)");
vector<double>::iterator iter = rArray.begin() + nIndex;
::std::nth_element( rArray.begin(), iter, rArray.end());
if (fDiff <= 0.0)
{
// Note: neg fDiff seen with forum-mso-en4-719754.xlsx with
// fPercentile of near 1 where approxFloor gave nIndex of nSize-1
// resulting in a non-zero tiny negative fDiff.
return *iter;
}
else
{
OSL_ENSURE(nIndex < nSize-1, "GetPercentile: wrong index(2)");
double fVal = *iter;
iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
return fVal + fDiff * (*iter - fVal);
}
}
}
double ScInterpreter::GetPercentileExclusive( vector<double> & rArray, double fPercentile )
{
size_t nSize1 = rArray.size() + 1;
if ( rArray.empty() || nSize1 == 1 || nGlobalError != FormulaError::NONE)
{
SetError( FormulaError::NoValue );
return 0.0;
}
if ( fPercentile * nSize1 < 1.0 || fPercentile * nSize1 > static_cast<double>( nSize1 - 1 ) )
{
SetError( FormulaError::IllegalParameter );
return 0.0;
}
size_t nIndex = static_cast<size_t>(::rtl::math::approxFloor( fPercentile * nSize1 - 1 ));
double fDiff = fPercentile * nSize1 - 1 - ::rtl::math::approxFloor( fPercentile * nSize1 - 1 );
OSL_ENSURE(nIndex < ( nSize1 - 1 ), "GetPercentile: wrong index(1)");
vector<double>::iterator iter = rArray.begin() + nIndex;
::std::nth_element( rArray.begin(), iter, rArray.end());
if (fDiff == 0.0)
return *iter;
else
{
OSL_ENSURE(nIndex < nSize1, "GetPercentile: wrong index(2)");
double fVal = *iter;
iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
return fVal + fDiff * (*iter - fVal);
}
}
void ScInterpreter::ScPercentile( bool bInclusive )
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double alpha = GetDouble();
if ( bInclusive ? ( alpha < 0.0 || alpha > 1.0 ) : ( alpha <= 0.0 || alpha >= 1.0 ) )
{
PushIllegalArgument();
return;
}
vector<double> aArray;
GetNumberSequenceArray( 1, aArray, false );
if ( aArray.empty() || nGlobalError != FormulaError::NONE )
{
PushNoValue();
return;
}
if ( bInclusive )
PushDouble( GetPercentile( aArray, alpha ));
else
PushDouble( GetPercentileExclusive( aArray, alpha ));
}
void ScInterpreter::ScQuartile( bool bInclusive )
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double fFlag = ::rtl::math::approxFloor(GetDouble());
if ( bInclusive ? ( fFlag < 0.0 || fFlag > 4.0 ) : ( fFlag <= 0.0 || fFlag >= 4.0 ) )
{
PushIllegalArgument();
return;
}
vector<double> aArray;
GetNumberSequenceArray( 1, aArray, false );
if ( aArray.empty() || nGlobalError != FormulaError::NONE )
{
PushNoValue();
return;
}
if ( bInclusive )
PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentile( aArray, 0.25 * fFlag ) );
else
PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentileExclusive( aArray, 0.25 * fFlag ) );
}
void ScInterpreter::ScModalValue()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCountMin( nParamCount, 1 ) )
return;
vector<double> aSortArray;
GetSortArray( nParamCount, aSortArray, nullptr, false, false );
SCSIZE nSize = aSortArray.size();
if (nSize == 0 || nGlobalError != FormulaError::NONE)
PushNoValue();
else
{
SCSIZE nMaxIndex = 0, nMax = 1, nCount = 1;
double nOldVal = aSortArray[0];
SCSIZE i;
for ( i = 1; i < nSize; i++)
{
if (aSortArray[i] == nOldVal)
nCount++;
else
{
nOldVal = aSortArray[i];
if (nCount > nMax)
{
nMax = nCount;
nMaxIndex = i-1;
}
nCount = 1;
}
}
if (nCount > nMax)
{
nMax = nCount;
nMaxIndex = i-1;
}
if (nMax == 1 && nCount == 1)
PushNoValue();
else if (nMax == 1)
PushDouble(nOldVal);
else
PushDouble(aSortArray[nMaxIndex]);
}
}
void ScInterpreter::ScModalValue_MS( bool bSingle )
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCountMin( nParamCount, 1 ) )
return;
vector<double> aArray;
GetNumberSequenceArray( nParamCount, aArray, false );
vector< double > aSortArray( aArray );
QuickSort( aSortArray, nullptr );
SCSIZE nSize = aSortArray.size();
if ( nSize == 0 || nGlobalError != FormulaError::NONE )
PushNoValue();
else
{
SCSIZE nMax = 1, nCount = 1;
double nOldVal = aSortArray[ 0 ];
vector< double > aResultArray( 1 );
SCSIZE i;
for ( i = 1; i < nSize; i++ )
{
if ( aSortArray[ i ] == nOldVal )
nCount++;
else
{
if ( nCount >= nMax && nCount > 1 )
{
if ( nCount > nMax )
{
nMax = nCount;
if ( aResultArray.size() != 1 )
vector< double >( 1 ).swap( aResultArray );
aResultArray[ 0 ] = nOldVal;
}
else
aResultArray.emplace_back( nOldVal );
}
nOldVal = aSortArray[ i ];
nCount = 1;
}
}
if ( nCount >= nMax && nCount > 1 )
{
if ( nCount > nMax )
vector< double >().swap( aResultArray );
aResultArray.emplace_back( nOldVal );
}
if ( nMax == 1 && nCount == 1 )
PushNoValue();
else if ( nMax == 1 )
PushDouble( nOldVal ); // there is only 1 result, no reordering needed
else
{
// sort resultArray according to ordering of aArray
vector< vector< double > > aOrder;
aOrder.resize( aResultArray.size(), vector< double >( 2 ) );
for ( i = 0; i < aResultArray.size(); i++ )
{
for ( SCSIZE j = 0; j < nSize; j++ )
{
if ( aArray[ j ] == aResultArray[ i ] )
{
aOrder[ i ][ 0 ] = aResultArray[ i ];
aOrder[ i ][ 1 ] = j;
break;
}
}
}
sort( aOrder.begin(), aOrder.end(), []( const std::vector< double >& lhs,
const std::vector< double >& rhs )
{ return lhs[ 1 ] < rhs[ 1 ]; } );
if ( bSingle )
PushDouble( aOrder[ 0 ][ 0 ] );
else
{
// put result in correct order in aResultArray
for ( i = 0; i < aResultArray.size(); i++ )
aResultArray[ i ] = aOrder[ i ][ 0 ];
ScMatrixRef pResMatrix = GetNewMat( 1, aResultArray.size(), true );
pResMatrix->PutDoubleVector( aResultArray, 0, 0 );
PushMatrix( pResMatrix );
}
}
}
}
void ScInterpreter::CalculateSmallLarge(bool bSmall)
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
SCSIZE nCol = 0, nRow = 0;
const auto aArray = GetRankNumberArray(nCol, nRow);
const size_t nRankArraySize = aArray.size();
if (nRankArraySize == 0 || nGlobalError != FormulaError::NONE)
{
PushNoValue();
return;
}
assert(nRankArraySize == nCol * nRow);
std::vector<SCSIZE> aRankArray;
aRankArray.reserve(nRankArraySize);
std::transform(aArray.begin(), aArray.end(), std::back_inserter(aRankArray),
[bSmall](double f) {
f = (bSmall ? rtl::math::approxFloor(f) : rtl::math::approxCeil(f));
// Valid ranks are >= 1.
if (f < 1.0 || !o3tl::convertsToAtMost(f, std::numeric_limits<SCSIZE>::max()))
return static_cast<SCSIZE>(0);
return static_cast<SCSIZE>(f);
});
vector<double> aSortArray;
GetNumberSequenceArray(1, aSortArray, false );
const SCSIZE nSize = aSortArray.size();
if (nSize == 0 || nGlobalError != FormulaError::NONE)
PushNoValue();
else if (nRankArraySize == 1)
{
const SCSIZE k = aRankArray[0];
if (k < 1 || nSize < k)
{
if (!std::isfinite(aArray[0]))
PushDouble(aArray[0]); // propagates error
else
PushNoValue();
}
else
{
vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k);
::std::nth_element( aSortArray.begin(), iPos, aSortArray.end());
PushDouble( *iPos);
}
}
else
{
std::set<SCSIZE> aIndices;
for (SCSIZE n : aRankArray)
{
if (1 <= n && n <= nSize)
aIndices.insert(bSmall ? n-1 : nSize-n);
}
// We can spare sorting when the total number of ranks is small enough.
// Find only the elements at given indices if, arbitrarily, the index size is
// smaller than 1/3 of the haystack array's size; just sort it squarely, otherwise.
if (aIndices.size() < nSize/3)
{
auto itBegin = aSortArray.begin();
for (SCSIZE i : aIndices)
{
auto it = aSortArray.begin() + i;
std::nth_element(itBegin, it, aSortArray.end());
itBegin = ++it;
}
}
else
std::sort(aSortArray.begin(), aSortArray.end());
std::vector<double> aResultArray;
aResultArray.reserve(nRankArraySize);
for (size_t i = 0; i < nRankArraySize; ++i)
{
const SCSIZE n = aRankArray[i];
if (1 <= n && n <= nSize)
aResultArray.push_back( aSortArray[bSmall ? n-1 : nSize-n]);
else if (!std::isfinite( aArray[i]))
aResultArray.push_back( aArray[i]); // propagate error
else
aResultArray.push_back( CreateDoubleError( FormulaError::IllegalArgument));
}
ScMatrixRef pResult = GetNewMat(nCol, nRow, aResultArray);
PushMatrix(pResult);
}
}
void ScInterpreter::ScLarge()
{
CalculateSmallLarge(false);
}
void ScInterpreter::ScSmall()
{
CalculateSmallLarge(true);
}
void ScInterpreter::ScPercentrank( bool bInclusive )
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
return;
double fSignificance = ( nParamCount == 3 ? ::rtl::math::approxFloor( GetDouble() ) : 3.0 );
if ( fSignificance < 1.0 )
{
PushIllegalArgument();
return;
}
double fNum = GetDouble();
vector<double> aSortArray;
GetSortArray( 1, aSortArray, nullptr, false, false );
SCSIZE nSize = aSortArray.size();
if ( nSize == 0 || nGlobalError != FormulaError::NONE )
PushNoValue();
else
{
if ( fNum < aSortArray[ 0 ] || fNum > aSortArray[ nSize - 1 ] )
PushNoValue();
else
{
double fRes;
if ( nSize == 1 )
fRes = 1.0; // fNum == aSortArray[ 0 ], see test above
else
fRes = GetPercentrank( aSortArray, fNum, bInclusive );
if ( fRes != 0.0 )
{
double fExp = ::rtl::math::approxFloor( log10( fRes ) ) + 1.0 - fSignificance;
fRes = ::rtl::math::round( fRes * pow( 10, -fExp ) ) / pow( 10, -fExp );
}
PushDouble( fRes );
}
}
}
double ScInterpreter::GetPercentrank( ::std::vector<double> & rArray, double fVal, bool bInclusive )
{
SCSIZE nSize = rArray.size();
double fRes;
if ( fVal == rArray[ 0 ] )
{
if ( bInclusive )
fRes = 0.0;
else
fRes = 1.0 / static_cast<double>( nSize + 1 );
}
else
{
SCSIZE nOldCount = 0;
double fOldVal = rArray[ 0 ];
SCSIZE i;
for ( i = 1; i < nSize && rArray[ i ] < fVal; i++ )
{
if ( rArray[ i ] != fOldVal )
{
nOldCount = i;
fOldVal = rArray[ i ];
}
}
if ( rArray[ i ] != fOldVal )
nOldCount = i;
if ( fVal == rArray[ i ] )
{
if ( bInclusive )
fRes = div( nOldCount, nSize - 1 );
else
fRes = static_cast<double>( i + 1 ) / static_cast<double>( nSize + 1 );
}
else
{
// nOldCount is the count of smaller entries
// fVal is between rArray[ nOldCount - 1 ] and rArray[ nOldCount ]
// use linear interpolation to find a position between the entries
if ( nOldCount == 0 )
{
OSL_FAIL( "should not happen" );
fRes = 0.0;
}
else
{
double fFract = ( fVal - rArray[ nOldCount - 1 ] ) /
( rArray[ nOldCount ] - rArray[ nOldCount - 1 ] );
if ( bInclusive )
fRes = div( static_cast<double>( nOldCount - 1 ) + fFract, nSize - 1 );
else
fRes = ( static_cast<double>(nOldCount) + fFract ) / static_cast<double>( nSize + 1 );
}
}
}
return fRes;
}
void ScInterpreter::ScTrimMean()
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double alpha = GetDouble();
if (alpha < 0.0 || alpha >= 1.0)
{
PushIllegalArgument();
return;
}
vector<double> aSortArray;
GetSortArray( 1, aSortArray, nullptr, false, false );
SCSIZE nSize = aSortArray.size();
if (nSize == 0 || nGlobalError != FormulaError::NONE)
PushNoValue();
else
{
sal_uLong nIndex = static_cast<sal_uLong>(::rtl::math::approxFloor(alpha*static_cast<double>(nSize)));
if (nIndex % 2 != 0)
nIndex--;
nIndex /= 2;
OSL_ENSURE(nIndex < nSize, "ScTrimMean: wrong index");
KahanSum fSum = 0.0;
for (SCSIZE i = nIndex; i < nSize-nIndex; i++)
fSum += aSortArray[i];
PushDouble(fSum.get()/static_cast<double>(nSize-2*nIndex));
}
}
std::vector<double> ScInterpreter::GetRankNumberArray( SCSIZE& rCol, SCSIZE& rRow )
{
std::vector<double> aArray;
switch (GetStackType())
{
case svDouble:
aArray.push_back(PopDouble());
rCol = rRow = 1;
break;
case svSingleRef:
{
ScAddress aAdr;
PopSingleRef(aAdr);
ScRefCellValue aCell(mrDoc, aAdr);
if (aCell.hasNumeric())
{
aArray.push_back(GetCellValue(aAdr, aCell));
rCol = rRow = 1;
}
}
break;
case svDoubleRef:
{
ScRange aRange;
PopDoubleRef(aRange, true);
if (nGlobalError != FormulaError::NONE)
break;
// give up unless the start and end are in the same sheet
if (aRange.aStart.Tab() != aRange.aEnd.Tab())
{
SetError(FormulaError::IllegalParameter);
break;
}
// the range already is in order
assert(aRange.aStart.Col() <= aRange.aEnd.Col());
assert(aRange.aStart.Row() <= aRange.aEnd.Row());
rCol = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
rRow = aRange.aEnd.Row() - aRange.aStart.Row() + 1;
aArray.reserve(rCol * rRow);
FormulaError nErr = FormulaError::NONE;
double fCellVal;
ScValueIterator aValIter(mrContext, aRange, mnSubTotalFlags);
if (aValIter.GetFirst(fCellVal, nErr))
{
do
aArray.push_back(fCellVal);
while (aValIter.GetNext(fCellVal, nErr) && nErr == FormulaError::NONE);
}
// Note that SMALL() and LARGE() rank parameters (2nd) have
// ParamClass::Value, so in array mode this is never hit and
// argument was converted to matrix instead, but for normal
// evaluation any non-numeric value including empty cell will
// result in error anyway, so just clear and propagate an existing
// error here already.
if (aArray.size() != rCol * rRow)
{
aArray.clear();
SetError(nErr);
}
}
break;
case svMatrix:
case svExternalSingleRef:
case svExternalDoubleRef:
{
ScMatrixRef pMat = GetMatrix();
if (!pMat)
break;
const SCSIZE nCount = pMat->GetElementCount();
aArray.reserve(nCount);
// Do not propagate errors from matrix elements as global error.
pMat->SetErrorInterpreter(nullptr);
if (pMat->IsNumeric())
{
for (SCSIZE i = 0; i < nCount; ++i)
aArray.push_back(pMat->GetDouble(i));
}
else
{
for (SCSIZE i = 0; i < nCount; ++i)
{
if (pMat->IsValue(i))
aArray.push_back( pMat->GetDouble(i));
else
aArray.push_back( CreateDoubleError( FormulaError::NoValue));
}
}
pMat->GetDimensions(rCol, rRow);
}
break;
default:
PopError();
SetError(FormulaError::IllegalParameter);
break;
}
return aArray;
}
void ScInterpreter::GetNumberSequenceArray( sal_uInt8 nParamCount, vector<double>& rArray, bool bConvertTextInArray )
{
ScAddress aAdr;
ScRange aRange;
const bool bIgnoreErrVal = bool(mnSubTotalFlags & SubtotalFlags::IgnoreErrVal);
short nParam = nParamCount;
size_t nRefInList = 0;
ReverseStack( nParamCount );
while (nParam-- > 0)
{
const StackVar eStackType = GetStackType();
switch (eStackType)
{
case svDouble :
rArray.push_back( PopDouble());
break;
case svSingleRef :
{
PopSingleRef( aAdr );
ScRefCellValue aCell(mrDoc, aAdr);
if (bIgnoreErrVal && aCell.hasError())
; // nothing
else if (aCell.hasNumeric())
rArray.push_back(GetCellValue(aAdr, aCell));
}
break;
case svDoubleRef :
case svRefList :
{
PopDoubleRef( aRange, nParam, nRefInList);
if (nGlobalError != FormulaError::NONE)
break;
aRange.PutInOrder();
SCSIZE nCellCount = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
nCellCount *= aRange.aEnd.Row() - aRange.aStart.Row() + 1;
rArray.reserve( rArray.size() + nCellCount);
FormulaError nErr = FormulaError::NONE;
double fCellVal;
ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
if (aValIter.GetFirst( fCellVal, nErr))
{
if (bIgnoreErrVal)
{
if (nErr == FormulaError::NONE)
rArray.push_back( fCellVal);
while (aValIter.GetNext( fCellVal, nErr))
{
if (nErr == FormulaError::NONE)
rArray.push_back( fCellVal);
}
}
else
{
rArray.push_back( fCellVal);
SetError(nErr);
while ((nErr == FormulaError::NONE) && aValIter.GetNext( fCellVal, nErr))
rArray.push_back( fCellVal);
SetError(nErr);
}
}
}
break;
case svMatrix :
case svExternalSingleRef:
case svExternalDoubleRef:
{
ScMatrixRef pMat = GetMatrix();
if (!pMat)
break;
SCSIZE nCount = pMat->GetElementCount();
rArray.reserve( rArray.size() + nCount);
if (pMat->IsNumeric())
{
if (bIgnoreErrVal)
{
for (SCSIZE i = 0; i < nCount; ++i)
{
const double fVal = pMat->GetDouble(i);
if (nGlobalError == FormulaError::NONE)
rArray.push_back( fVal);
else
nGlobalError = FormulaError::NONE;
}
}
else
{
for (SCSIZE i = 0; i < nCount; ++i)
rArray.push_back( pMat->GetDouble(i));
}
}
else if (bConvertTextInArray && eStackType == svMatrix)
{
for (SCSIZE i = 0; i < nCount; ++i)
{
if ( pMat->IsValue( i ) )
{
if (bIgnoreErrVal)
{
const double fVal = pMat->GetDouble(i);
if (nGlobalError == FormulaError::NONE)
rArray.push_back( fVal);
else
nGlobalError = FormulaError::NONE;
}
else
rArray.push_back( pMat->GetDouble(i));
}
else
{
// tdf#88547 try to convert string to (date)value
OUString aStr = pMat->GetString( i ).getString();
if ( aStr.getLength() > 0 )
{
FormulaError nErr = nGlobalError;
nGlobalError = FormulaError::NONE;
double fVal = ConvertStringToValue( aStr );
if ( nGlobalError == FormulaError::NONE )
{
rArray.push_back( fVal );
nGlobalError = nErr;
}
else
{
if (!bIgnoreErrVal)
rArray.push_back( CreateDoubleError( FormulaError::NoValue));
// Propagate previous error if any, else
// the current #VALUE! error, unless
// ignoring error values.
if (nErr != FormulaError::NONE)
nGlobalError = nErr;
else if (!bIgnoreErrVal)
nGlobalError = FormulaError::NoValue;
else
nGlobalError = FormulaError::NONE;
}
}
}
}
}
else
{
if (bIgnoreErrVal)
{
for (SCSIZE i = 0; i < nCount; ++i)
{
if (pMat->IsValue(i))
{
const double fVal = pMat->GetDouble(i);
if (nGlobalError == FormulaError::NONE)
rArray.push_back( fVal);
else
nGlobalError = FormulaError::NONE;
}
}
}
else
{
for (SCSIZE i = 0; i < nCount; ++i)
{
if (pMat->IsValue(i))
rArray.push_back( pMat->GetDouble(i));
}
}
}
}
break;
default :
PopError();
SetError( FormulaError::IllegalParameter);
break;
}
if (nGlobalError != FormulaError::NONE)
break; // while
}
// nParam > 0 in case of error, clean stack environment and obtain earlier
// error if there was one.
while (nParam-- > 0)
PopError();
}
void ScInterpreter::DecoladeRow( ScSortInfoArray* pArray, SCROW nRow1, SCROW nRow2 )
{
SCROW nRow;
int nMax = nRow2 - nRow1;
for (SCROW i = nRow1; (i + 4) <= nRow2; i += 4)
{
nRow = comphelper::rng::uniform_int_distribution(0, nMax - 1);
pArray->Swap(i, nRow1 + nRow);
}
}
std::unique_ptr<ScSortInfoArray> ScInterpreter::CreateFastSortInfoArray(
const ScSortParam& rSortParam, bool bMatrix, SCCOLROW nInd1, SCCOLROW nInd2 )
{
sal_uInt16 nUsedSorts = 1;
while (nUsedSorts < rSortParam.GetSortKeyCount() && rSortParam.maKeyState[nUsedSorts].bDoSort)
nUsedSorts++;
std::unique_ptr<ScSortInfoArray> pArray(new ScSortInfoArray(nUsedSorts, nInd1, nInd2));
if (rSortParam.bByRow)
{
for (sal_uInt16 nSort = 0; nSort < nUsedSorts; nSort++)
{
if (!bMatrix)
{
SCCOL nCol = static_cast<SCCOL>(rSortParam.maKeyState[nSort].nField);
std::optional<sc::ColumnIterator> pIter = mrDoc.GetColumnIterator(rSortParam.nSourceTab, nCol, nInd1, nInd2);
assert(pIter->hasCell());
for (SCROW nRow = nInd1; nRow <= nInd2; nRow++, pIter->next())
{
ScSortInfo& rInfo = pArray->Get(nSort, nRow);
rInfo.maCell = pIter->getCell();
rInfo.nOrg = nRow;
}
}
else
{
for (SCROW nRow = nInd1; nRow <= nInd2; nRow++)
{
ScSortInfo& rInfo = pArray->Get(nSort, nRow);
rInfo.nOrg = nRow;
}
}
}
}
else
{
for (sal_uInt16 nSort = 0; nSort < nUsedSorts; nSort++)
{
if (!bMatrix)
{
SCROW nRow = rSortParam.maKeyState[nSort].nField;
for (SCCOL nCol = static_cast<SCCOL>(nInd1);
nCol <= static_cast<SCCOL>(nInd2); nCol++)
{
ScSortInfo& rInfo = pArray->Get(nSort, nCol);
rInfo.maCell = mrDoc.GetRefCellValue(ScAddress(nCol, nRow, rSortParam.nSourceTab));
rInfo.nOrg = nCol;
}
}
else
{
for (SCCOL nCol = static_cast<SCCOL>(nInd1);
nCol <= static_cast<SCCOL>(nInd2); nCol++)
{
ScSortInfo& rInfo = pArray->Get(nSort, nCol);
rInfo.nOrg = nCol;
}
}
}
}
return pArray;
}
std::vector<SCCOLROW> ScInterpreter::GetSortOrder( const ScSortParam& rSortParam, const ScMatrixRef& pMatSrc )
{
std::vector<SCCOLROW> aOrderIndices;
aSortParam = rSortParam;
if (rSortParam.bByRow)
{
const SCROW nLastRow = rSortParam.nRow2;
const SCROW nRow1 = (rSortParam.bHasHeader ? rSortParam.nRow1 + 1 : rSortParam.nRow1);
if (nRow1 < nLastRow)
{
std::unique_ptr<ScSortInfoArray> pArray(CreateFastSortInfoArray(
aSortParam, (pMatSrc != nullptr), nRow1, nLastRow));
if (nLastRow - nRow1 > 255)
DecoladeRow(pArray.get(), nRow1, nLastRow);
QuickSort(pArray.get(), pMatSrc, nRow1, nLastRow);
aOrderIndices = pArray->GetOrderIndices();
}
}
else
{
const SCCOL nLastCol = rSortParam.nCol2;
const SCCOL nCol1 = (rSortParam.bHasHeader ? rSortParam.nCol1 + 1 : rSortParam.nCol1);
if (nCol1 < nLastCol)
{
std::unique_ptr<ScSortInfoArray> pArray(CreateFastSortInfoArray(
aSortParam, (pMatSrc != nullptr), nCol1, nLastCol));
QuickSort(pArray.get(), pMatSrc, nCol1, nLastCol);
aOrderIndices = pArray->GetOrderIndices();
}
}
return aOrderIndices;
}
ScMatrixRef ScInterpreter::CreateSortedMatrix( const ScSortParam& rSortParam, const ScMatrixRef& pMatSrc,
const ScRange& rSourceRange, const std::vector<SCCOLROW>& rSortArray, SCSIZE nsC, SCSIZE nsR )
{
SCCOLROW nStartPos = (!rSortParam.bByRow ? rSortParam.nCol1 : rSortParam.nRow1);
size_t nCount = rSortArray.size();
std::vector<SCCOLROW> aPosTable(nCount);
for (size_t i = 0; i < nCount; ++i)
aPosTable[rSortArray[i] - nStartPos] = i;
ScMatrixRef pResMat = nullptr;
if (!rSortArray.empty())
{
pResMat = GetNewMat(nsC, nsR, /*bEmpty*/true);
if (!pMatSrc)
{
ScCellIterator aCellIter(mrDoc, rSourceRange);
for (bool bHas = aCellIter.first(); bHas; bHas = aCellIter.next())
{
SCSIZE nThisCol = static_cast<SCSIZE>(aCellIter.GetPos().Col() - rSourceRange.aStart.Col());
SCSIZE nThisRow = static_cast<SCSIZE>(aCellIter.GetPos().Row() - rSourceRange.aStart.Row());
ScRefCellValue aCell = aCellIter.getRefCellValue();
if (aCell.hasNumeric())
{
if (rSortParam.bByRow)
pResMat->PutDouble(GetCellValue(aCellIter.GetPos(), aCell), nThisCol, aPosTable[nThisRow]);
else
pResMat->PutDouble(GetCellValue(aCellIter.GetPos(), aCell), aPosTable[nThisCol], nThisRow);
}
else
{
svl::SharedString aStr;
GetCellString(aStr, aCell);
if (rSortParam.bByRow)
pResMat->PutString(aStr, nThisCol, aPosTable[nThisRow]);
else
pResMat->PutString(aStr, aPosTable[nThisCol], nThisRow);
}
}
}
else
{
for (SCCOL ci = rSourceRange.aStart.Col(); ci <= rSourceRange.aEnd.Col(); ci++)
{
for (SCROW rj = rSourceRange.aStart.Row(); rj <= rSourceRange.aEnd.Row(); rj++)
{
if (pMatSrc->IsEmptyCell(ci, rj))
{
if (rSortParam.bByRow)
pResMat->PutEmpty(ci, aPosTable[rj]);
else
pResMat->PutEmpty(aPosTable[ci], rj);
}
else if (pMatSrc->IsStringOrEmpty(ci, rj))
{
if (rSortParam.bByRow)
pResMat->PutString(pMatSrc->GetString(ci, rj), ci, aPosTable[rj]);
else
pResMat->PutString(pMatSrc->GetString(ci, rj), aPosTable[ci], rj);
}
else
{
if (rSortParam.bByRow)
pResMat->PutDouble(pMatSrc->GetDouble(ci, rj), ci, aPosTable[rj]);
else
pResMat->PutDouble(pMatSrc->GetDouble(ci, rj), aPosTable[ci], rj);
}
}
}
}
}
return pResMat;
}
void ScInterpreter::QuickSort( ScSortInfoArray* pArray, const ScMatrixRef& pMatSrc, SCCOLROW nLo, SCCOLROW nHi )
{
if ((nHi - nLo) == 1)
{
if (Compare(pArray, pMatSrc, nLo, nHi) > 0)
pArray->Swap( nLo, nHi );
}
else
{
SCCOLROW ni = nLo;
SCCOLROW nj = nHi;
do
{
while ((ni <= nHi) && (Compare(pArray, pMatSrc, ni, nLo)) < 0)
ni++;
while ((nj >= nLo) && (Compare(pArray, pMatSrc, nLo, nj)) < 0)
nj--;
if (ni <= nj)
{
if (ni != nj)
pArray->Swap( ni, nj );
ni++;
nj--;
}
} while (ni < nj);
if ((nj - nLo) < (nHi - ni))
{
if (nLo < nj)
QuickSort(pArray, pMatSrc, nLo, nj);
if (ni < nHi)
QuickSort(pArray, pMatSrc, ni, nHi);
}
else
{
if (ni < nHi)
QuickSort(pArray, pMatSrc, ni, nHi);
if (nLo < nj)
QuickSort(pArray, pMatSrc, nLo, nj);
}
}
}
short ScInterpreter::Compare( ScSortInfoArray* pArray, const ScMatrixRef& pMatSrc, SCCOLROW nIndex1, SCCOLROW nIndex2 ) const
{
short nRes;
sal_uInt16 nSort = 0;
do
{
ScSortInfo& rInfo1 = pArray->Get( nSort, nIndex1 );
ScSortInfo& rInfo2 = pArray->Get( nSort, nIndex2 );
if (!pMatSrc)
{
nRes = CompareCell(nSort, rInfo1.maCell, rInfo2.maCell);
}
else
{
if (aSortParam.bByRow)
nRes = CompareMatrixCell( pMatSrc, nSort,
static_cast<SCCOL>(aSortParam.maKeyState[nSort].nField), rInfo1.nOrg,
static_cast<SCCOL>(aSortParam.maKeyState[nSort].nField), rInfo2.nOrg );
else
nRes = CompareMatrixCell( pMatSrc, nSort,
static_cast<SCCOL>(rInfo1.nOrg), aSortParam.maKeyState[nSort].nField,
static_cast<SCCOL>(rInfo2.nOrg), aSortParam.maKeyState[nSort].nField );
}
} while ( nRes == 0 && ++nSort < pArray->GetUsedSorts() );
if( nRes == 0 )
{
ScSortInfo& rInfo1 = pArray->Get( 0, nIndex1 );
ScSortInfo& rInfo2 = pArray->Get( 0, nIndex2 );
if( rInfo1.nOrg < rInfo2.nOrg )
nRes = -1;
else if( rInfo1.nOrg > rInfo2.nOrg )
nRes = 1;
}
return nRes;
}
short ScInterpreter::CompareCell( sal_uInt16 nSort,
ScRefCellValue& rCell1, ScRefCellValue& rCell2 ) const
{
short nRes = 0;
CellType eType1 = rCell1.getType(), eType2 = rCell2.getType();
if (!rCell1.isEmpty())
{
if (!rCell2.isEmpty())
{
bool bErr1 = false;
bool bStr1 = ( eType1 != CELLTYPE_VALUE );
if (eType1 == CELLTYPE_FORMULA)
{
if (rCell1.getFormula()->GetErrCode() != FormulaError::NONE)
{
bErr1 = true;
bStr1 = false;
}
else if (rCell1.getFormula()->IsValue())
{
bStr1 = false;
}
}
bool bErr2 = false;
bool bStr2 = ( eType2 != CELLTYPE_VALUE );
if (eType2 == CELLTYPE_FORMULA)
{
if (rCell2.getFormula()->GetErrCode() != FormulaError::NONE)
{
bErr2 = true;
bStr2 = false;
}
else if (rCell2.getFormula()->IsValue())
{
bStr2 = false;
}
}
if ( bStr1 && bStr2 ) // only compare strings as strings!
{
OUString aStr1;
OUString aStr2;
if (eType1 == CELLTYPE_STRING)
aStr1 = rCell1.getSharedString()->getString();
else
aStr1 = rCell1.getString(&mrDoc);
if (eType2 == CELLTYPE_STRING)
aStr2 = rCell2.getSharedString()->getString();
else
aStr2 = rCell2.getString(&mrDoc);
CollatorWrapper& rSortCollator = ScGlobal::GetCollator(aSortParam.bCaseSens);
nRes = static_cast<short>( rSortCollator.compareString( aStr1, aStr2 ) );
}
else if ( bStr1 ) // String <-> Number or Error
{
if (bErr2)
nRes = -1; // String in front of Error
else
nRes = 1; // Number in front of String
}
else if ( bStr2 ) // Number or Error <-> String
{
if (bErr1)
nRes = 1; // String in front of Error
else
nRes = -1; // Number in front of String
}
else if (bErr1 && bErr2)
{
// nothing, two Errors are equal
}
else if (bErr1) // Error <-> Number
{
nRes = 1; // Number in front of Error
}
else if (bErr2) // Number <-> Error
{
nRes = -1; // Number in front of Error
}
else // Mixed numbers
{
double nVal1 = rCell1.getValue();
double nVal2 = rCell2.getValue();
if (nVal1 < nVal2)
nRes = -1;
else if (nVal1 > nVal2)
nRes = 1;
}
if ( !aSortParam.maKeyState[nSort].bAscending )
nRes = -nRes;
}
else
nRes = -1;
}
else
{
if (!rCell2.isEmpty())
nRes = 1;
else
nRes = 0; // both empty
}
return nRes;
}
short ScInterpreter::CompareMatrixCell( const ScMatrixRef& pMatSrc, sal_uInt16 nSort, SCCOL nCell1Col, SCROW nCell1Row,
SCCOL nCell2Col, SCROW nCell2Row ) const
{
short nRes = 0;
// 1st value
bool bCell1Empty = false;
bool bCell1Value = false;
if (pMatSrc->IsEmpty(nCell1Col, nCell1Row))
bCell1Empty = true;
else if (pMatSrc->IsStringOrEmpty(nCell1Col, nCell1Row))
bCell1Value = false;
else
bCell1Value = true;
// 2nd value
bool bCell2Empty = false;
bool bCell2Value = false;
if (pMatSrc->IsEmpty(nCell2Col, nCell2Row))
bCell2Empty = true;
else if (pMatSrc->IsStringOrEmpty(nCell2Col, nCell2Row))
bCell2Value = false;
else
bCell2Value = true;
if (!bCell1Empty)
{
if (!bCell2Empty)
{
if ( !bCell1Value && !bCell2Value ) // only compare strings as strings!
{
OUString aStr1 = pMatSrc->GetString(nCell1Col, nCell1Row).getString();
OUString aStr2 = pMatSrc->GetString(nCell2Col, nCell2Row).getString();
CollatorWrapper& rSortCollator = ScGlobal::GetCollator(aSortParam.bCaseSens);
nRes = static_cast<short>( rSortCollator.compareString( aStr1, aStr2 ) );
}
else if ( !bCell1Value ) // String <-> Number or Error
{
nRes = 1; // Number in front of String
}
else if ( !bCell2Value ) // Number or Error <-> String
{
nRes = -1; // Number in front of String
}
else // Mixed numbers
{
double nVal1 = pMatSrc->GetDouble(nCell1Col, nCell1Row);
double nVal2 = pMatSrc->GetDouble(nCell2Col, nCell2Row);
if (nVal1 < nVal2)
nRes = -1;
else if (nVal1 > nVal2)
nRes = 1;
}
if ( !aSortParam.maKeyState[nSort].bAscending )
nRes = -nRes;
}
else
nRes = -1;
}
else
{
if (!bCell2Empty)
nRes = 1;
else
nRes = 0; // both empty
}
return nRes;
}
void ScInterpreter::GetSortArray( sal_uInt8 nParamCount, vector<double>& rSortArray, vector<tools::Long>* pIndexOrder, bool bConvertTextInArray, bool bAllowEmptyArray )
{
GetNumberSequenceArray( nParamCount, rSortArray, bConvertTextInArray );
if (rSortArray.size() > MAX_COUNT_DOUBLE_FOR_SORT(mrDoc.GetSheetLimits()))
SetError( FormulaError::MatrixSize);
else if ( rSortArray.empty() )
{
if ( bAllowEmptyArray )
return;
SetError( FormulaError::NoValue);
}
if (nGlobalError == FormulaError::NONE)
QuickSort( rSortArray, pIndexOrder);
}
static void lcl_QuickSort( tools::Long nLo, tools::Long nHi, vector<double>& rSortArray, vector<tools::Long>* pIndexOrder )
{
// If pIndexOrder is not NULL, we assume rSortArray.size() == pIndexOrder->size().
using ::std::swap;
if (nHi - nLo == 1)
{
if (rSortArray[nLo] > rSortArray[nHi])
{
swap(rSortArray[nLo], rSortArray[nHi]);
if (pIndexOrder)
swap(pIndexOrder->at(nLo), pIndexOrder->at(nHi));
}
return;
}
tools::Long ni = nLo;
tools::Long nj = nHi;
do
{
double fLo = rSortArray[nLo];
while (ni <= nHi && rSortArray[ni] < fLo) ni++;
while (nj >= nLo && fLo < rSortArray[nj]) nj--;
if (ni <= nj)
{
if (ni != nj)
{
swap(rSortArray[ni], rSortArray[nj]);
if (pIndexOrder)
swap(pIndexOrder->at(ni), pIndexOrder->at(nj));
}
++ni;
--nj;
}
}
while (ni < nj);
if ((nj - nLo) < (nHi - ni))
{
if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
}
else
{
if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
}
}
void ScInterpreter::QuickSort( vector<double>& rSortArray, vector<tools::Long>* pIndexOrder )
{
tools::Long n = static_cast<tools::Long>(rSortArray.size());
if (pIndexOrder)
{
pIndexOrder->clear();
pIndexOrder->reserve(n);
for (tools::Long i = 0; i < n; ++i)
pIndexOrder->push_back(i);
}
if (n < 2)
return;
size_t nValCount = rSortArray.size();
for (size_t i = 0; (i + 4) <= nValCount-1; i += 4)
{
size_t nInd = comphelper::rng::uniform_size_distribution(0, nValCount-2);
::std::swap( rSortArray[i], rSortArray[nInd]);
if (pIndexOrder)
::std::swap( pIndexOrder->at(i), pIndexOrder->at(nInd));
}
lcl_QuickSort(0, n-1, rSortArray, pIndexOrder);
}
void ScInterpreter::ScRank( bool bAverage )
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
return;
bool bAscending;
if ( nParamCount == 3 )
bAscending = GetBool();
else
bAscending = false;
vector<double> aSortArray;
GetSortArray( 1, aSortArray, nullptr, false, false );
double fVal = GetDouble();
SCSIZE nSize = aSortArray.size();
if ( nSize == 0 || nGlobalError != FormulaError::NONE )
PushNoValue();
else
{
if ( fVal < aSortArray[ 0 ] || fVal > aSortArray[ nSize - 1 ] )
PushError( FormulaError::NotAvailable);
else
{
double fLastPos = 0;
double fFirstPos = -1.0;
bool bFinished = false;
SCSIZE i;
for (i = 0; i < nSize && !bFinished; i++)
{
if ( aSortArray[ i ] == fVal )
{
if ( fFirstPos < 0 )
fFirstPos = i + 1.0;
}
else
{
if ( aSortArray[ i ] > fVal )
{
fLastPos = i;
bFinished = true;
}
}
}
if ( !bFinished )
fLastPos = i;
if (fFirstPos <= 0)
{
PushError( FormulaError::NotAvailable);
}
else if ( !bAverage )
{
if ( bAscending )
PushDouble( fFirstPos );
else
PushDouble( nSize + 1.0 - fLastPos );
}
else
{
if ( bAscending )
PushDouble( ( fFirstPos + fLastPos ) / 2.0 );
else
PushDouble( nSize + 1.0 - ( fFirstPos + fLastPos ) / 2.0 );
}
}
}
}
void ScInterpreter::ScAveDev()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCountMin( nParamCount, 1 ) )
return;
sal_uInt16 SaveSP = sp;
double nMiddle = 0.0;
KahanSum rVal = 0.0;
double rValCount = 0.0;
ScAddress aAdr;
ScRange aRange;
short nParam = nParamCount;
size_t nRefInList = 0;
while (nParam-- > 0)
{
switch (GetStackType())
{
case svDouble :
rVal += GetDouble();
rValCount++;
break;
case svSingleRef :
{
PopSingleRef( aAdr );
ScRefCellValue aCell(mrDoc, aAdr);
if (aCell.hasNumeric())
{
rVal += GetCellValue(aAdr, aCell);
rValCount++;
}
}
break;
case svDoubleRef :
case svRefList :
{
FormulaError nErr = FormulaError::NONE;
double nCellVal;
PopDoubleRef( aRange, nParam, nRefInList);
ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
if (aValIter.GetFirst(nCellVal, nErr))
{
rVal += nCellVal;
rValCount++;
SetError(nErr);
while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
{
rVal += nCellVal;
rValCount++;
}
SetError(nErr);
}
}
break;
case svMatrix :
case svExternalSingleRef:
case svExternalDoubleRef:
{
ScMatrixRef pMat = GetMatrix();
if (pMat)
{
SCSIZE nCount = pMat->GetElementCount();
if (pMat->IsNumeric())
{
for (SCSIZE nElem = 0; nElem < nCount; nElem++)
{
rVal += pMat->GetDouble(nElem);
rValCount++;
}
}
else
{
for (SCSIZE nElem = 0; nElem < nCount; nElem++)
if (!pMat->IsStringOrEmpty(nElem))
{
rVal += pMat->GetDouble(nElem);
rValCount++;
}
}
}
}
break;
default :
SetError(FormulaError::IllegalParameter);
break;
}
}
if (nGlobalError != FormulaError::NONE)
{
PushError( nGlobalError);
return;
}
nMiddle = rVal.get() / rValCount;
sp = SaveSP;
rVal = 0.0;
nParam = nParamCount;
nRefInList = 0;
while (nParam-- > 0)
{
switch (GetStackType())
{
case svDouble :
rVal += std::abs(GetDouble() - nMiddle);
break;
case svSingleRef :
{
PopSingleRef( aAdr );
ScRefCellValue aCell(mrDoc, aAdr);
if (aCell.hasNumeric())
rVal += std::abs(GetCellValue(aAdr, aCell) - nMiddle);
}
break;
case svDoubleRef :
case svRefList :
{
FormulaError nErr = FormulaError::NONE;
double nCellVal;
PopDoubleRef( aRange, nParam, nRefInList);
ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
if (aValIter.GetFirst(nCellVal, nErr))
{
rVal += std::abs(nCellVal - nMiddle);
while (aValIter.GetNext(nCellVal, nErr))
rVal += std::abs(nCellVal - nMiddle);
}
}
break;
case svMatrix :
case svExternalSingleRef:
case svExternalDoubleRef:
{
ScMatrixRef pMat = GetMatrix();
if (pMat)
{
SCSIZE nCount = pMat->GetElementCount();
if (pMat->IsNumeric())
{
for (SCSIZE nElem = 0; nElem < nCount; nElem++)
{
rVal += std::abs(pMat->GetDouble(nElem) - nMiddle);
}
}
else
{
for (SCSIZE nElem = 0; nElem < nCount; nElem++)
{
if (!pMat->IsStringOrEmpty(nElem))
rVal += std::abs(pMat->GetDouble(nElem) - nMiddle);
}
}
}
}
break;
default : SetError(FormulaError::IllegalParameter); break;
}
}
PushDouble(rVal.get() / rValCount);
}
void ScInterpreter::ScDevSq()
{
auto VarResult = []( double fVal, size_t /*nValCount*/ )
{
return fVal;
};
GetStVarParams( false /*bTextAsZero*/, VarResult);
}
void ScInterpreter::ScProbability()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
return;
double fUp, fLo;
fUp = GetDouble();
if (nParamCount == 4)
fLo = GetDouble();
else
fLo = fUp;
if (fLo > fUp)
std::swap( fLo, fUp );
ScMatrixRef pMatP = GetMatrix();
ScMatrixRef pMatW = GetMatrix();
if (!pMatP || !pMatW)
PushIllegalParameter();
else
{
SCSIZE nC1, nC2;
SCSIZE nR1, nR2;
pMatP->GetDimensions(nC1, nR1);
pMatW->GetDimensions(nC2, nR2);
if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 ||
nC2 == 0 || nR2 == 0)
PushNA();
else
{
KahanSum fSum = 0.0;
KahanSum fRes = 0.0;
bool bStop = false;
double fP, fW;
for ( SCSIZE i = 0; i < nC1 && !bStop; i++ )
{
for (SCSIZE j = 0; j < nR1 && !bStop; ++j )
{
if (pMatP->IsValue(i,j) && pMatW->IsValue(i,j))
{
fP = pMatP->GetDouble(i,j);
fW = pMatW->GetDouble(i,j);
if (fP < 0.0 || fP > 1.0)
bStop = true;
else
{
fSum += fP;
if (fW >= fLo && fW <= fUp)
fRes += fP;
}
}
else
SetError( FormulaError::IllegalArgument);
}
}
if (bStop || std::abs((fSum -1.0).get()) > 1.0E-7)
PushNoValue();
else
PushDouble(fRes.get());
}
}
}
void ScInterpreter::ScCorrel()
{
// This is identical to ScPearson()
ScPearson();
}
void ScInterpreter::ScCovarianceP()
{
CalculatePearsonCovar( false, false, false );
}
void ScInterpreter::ScCovarianceS()
{
CalculatePearsonCovar( false, false, true );
}
void ScInterpreter::ScPearson()
{
CalculatePearsonCovar( true, false, false );
}
void ScInterpreter::CalculatePearsonCovar( bool _bPearson, bool _bStexy, bool _bSample )
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
ScMatrixRef pMat1 = GetMatrix();
ScMatrixRef pMat2 = GetMatrix();
if (!pMat1 || !pMat2)
{
PushIllegalParameter();
return;
}
SCSIZE nC1, nC2;
SCSIZE nR1, nR2;
pMat1->GetDimensions(nC1, nR1);
pMat2->GetDimensions(nC2, nR2);
if (nR1 != nR2 || nC1 != nC2)
{
PushIllegalArgument();
return;
}
/* #i78250#
* (sum((X-MeanX)(Y-MeanY)))/N equals (SumXY)/N-MeanX*MeanY mathematically,
* but the latter produces wrong results if the absolute values are high,
* for example above 10^8
*/
double fCount = 0.0;
KahanSum fSumX = 0.0;
KahanSum fSumY = 0.0;
for (SCSIZE i = 0; i < nC1; i++)
{
for (SCSIZE j = 0; j < nR1; j++)
{
if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
{
fSumX += pMat1->GetDouble(i,j);
fSumY += pMat2->GetDouble(i,j);
fCount++;
}
}
}
if (fCount < (_bStexy ? 3.0 : (_bSample ? 2.0 : 1.0)))
PushNoValue();
else
{
KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
KahanSum fSumSqrDeltaY = 0.0; // sum of (ValY-MeanY)^2
const double fMeanX = fSumX.get() / fCount;
const double fMeanY = fSumY.get() / fCount;
for (SCSIZE i = 0; i < nC1; i++)
{
for (SCSIZE j = 0; j < nR1; j++)
{
if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
{
const double fValX = pMat1->GetDouble(i,j);
const double fValY = pMat2->GetDouble(i,j);
fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
if ( _bPearson )
{
fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
fSumSqrDeltaY += (fValY - fMeanY) * (fValY - fMeanY);
}
}
}
}
if ( _bPearson )
{
// tdf#94962 - Values below the numerical limit lead to unexpected results
if (fSumSqrDeltaX < ::std::numeric_limits<double>::min()
|| (!_bStexy && fSumSqrDeltaY < ::std::numeric_limits<double>::min()))
PushError( FormulaError::DivisionByZero);
else if ( _bStexy )
PushDouble( sqrt( ( fSumSqrDeltaY - fSumDeltaXDeltaY *
fSumDeltaXDeltaY / fSumSqrDeltaX ).get() / (fCount-2)));
else
PushDouble( fSumDeltaXDeltaY.get() / sqrt( fSumSqrDeltaX.get() * fSumSqrDeltaY.get() ));
}
else
{
if ( _bSample )
PushDouble( fSumDeltaXDeltaY.get() / ( fCount - 1 ) );
else
PushDouble( fSumDeltaXDeltaY.get() / fCount);
}
}
}
void ScInterpreter::ScRSQ()
{
// Same as ScPearson()*ScPearson()
ScPearson();
if (nGlobalError != FormulaError::NONE)
return;
switch (GetStackType())
{
case svDouble:
{
double fVal = PopDouble();
PushDouble( fVal * fVal);
}
break;
default:
PopError();
PushNoValue();
}
}
void ScInterpreter::ScSTEYX()
{
CalculatePearsonCovar( true, true, false );
}
void ScInterpreter::CalculateSlopeIntercept(bool bSlope)
{
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
ScMatrixRef pMat1 = GetMatrix();
ScMatrixRef pMat2 = GetMatrix();
if (!pMat1 || !pMat2)
{
PushIllegalParameter();
return;
}
SCSIZE nC1, nC2;
SCSIZE nR1, nR2;
pMat1->GetDimensions(nC1, nR1);
pMat2->GetDimensions(nC2, nR2);
if (nR1 != nR2 || nC1 != nC2)
{
PushIllegalArgument();
return;
}
// #i78250# numerical stability improved
double fCount = 0.0;
KahanSum fSumX = 0.0;
KahanSum fSumY = 0.0;
for (SCSIZE i = 0; i < nC1; i++)
{
for (SCSIZE j = 0; j < nR1; j++)
{
if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
{
fSumX += pMat1->GetDouble(i,j);
fSumY += pMat2->GetDouble(i,j);
fCount++;
}
}
}
if (fCount < 1.0)
PushNoValue();
else
{
KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
double fMeanX = fSumX.get() / fCount;
double fMeanY = fSumY.get() / fCount;
for (SCSIZE i = 0; i < nC1; i++)
{
for (SCSIZE j = 0; j < nR1; j++)
{
if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
{
double fValX = pMat1->GetDouble(i,j);
double fValY = pMat2->GetDouble(i,j);
fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
}
}
}
if (fSumSqrDeltaX == 0.0)
PushError( FormulaError::DivisionByZero);
else
{
if ( bSlope )
PushDouble( fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get());
else
PushDouble( fMeanY - fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get() * fMeanX);
}
}
}
void ScInterpreter::ScSlope()
{
CalculateSlopeIntercept(true);
}
void ScInterpreter::ScIntercept()
{
CalculateSlopeIntercept(false);
}
void ScInterpreter::ScForecast()
{
if ( !MustHaveParamCount( GetByte(), 3 ) )
return;
ScMatrixRef pMat1 = GetMatrix();
ScMatrixRef pMat2 = GetMatrix();
if (!pMat1 || !pMat2)
{
PushIllegalParameter();
return;
}
SCSIZE nC1, nC2;
SCSIZE nR1, nR2;
pMat1->GetDimensions(nC1, nR1);
pMat2->GetDimensions(nC2, nR2);
if (nR1 != nR2 || nC1 != nC2)
{
PushIllegalArgument();
return;
}
double fVal = GetDouble();
// #i78250# numerical stability improved
double fCount = 0.0;
KahanSum fSumX = 0.0;
KahanSum fSumY = 0.0;
for (SCSIZE i = 0; i < nC1; i++)
{
for (SCSIZE j = 0; j < nR1; j++)
{
if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
{
fSumX += pMat1->GetDouble(i,j);
fSumY += pMat2->GetDouble(i,j);
fCount++;
}
}
}
if (fCount < 1.0)
PushNoValue();
else
{
KahanSum fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
KahanSum fSumSqrDeltaX = 0.0; // sum of (ValX-MeanX)^2
double fMeanX = fSumX.get() / fCount;
double fMeanY = fSumY.get() / fCount;
for (SCSIZE i = 0; i < nC1; i++)
{
for (SCSIZE j = 0; j < nR1; j++)
{
if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
{
double fValX = pMat1->GetDouble(i,j);
double fValY = pMat2->GetDouble(i,j);
fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
fSumSqrDeltaX += (fValX - fMeanX) * (fValX - fMeanX);
}
}
}
if (fSumSqrDeltaX == 0.0)
PushError( FormulaError::DivisionByZero);
else
PushDouble( fMeanY + fSumDeltaXDeltaY.get() / fSumSqrDeltaX.get() * (fVal - fMeanX));
}
}
static void lcl_roundUpNearestPow2(SCSIZE& nNum, SCSIZE& nNumBits)
{
// Find the least power of 2 that is less than or equal to nNum.
SCSIZE nPow2(1);
nNumBits = std::numeric_limits<SCSIZE>::digits;
nPow2 <<= (nNumBits - 1);
while (nPow2 >= 1)
{
if (nNum & nPow2)
break;
--nNumBits;
nPow2 >>= 1;
}
if (nPow2 != nNum)
{
assert(nPow2 < 1UL << (std::numeric_limits<unsigned long>::digits - 1));
nNum = nPow2 ? (nPow2 << 1) : 1;
}
else
--nNumBits;
}
static SCSIZE lcl_bitReverse(SCSIZE nIn, SCSIZE nBound)
{
SCSIZE nOut = 0;
for (SCSIZE nMask = 1; nMask < nBound; nMask <<= 1)
{
nOut <<= 1;
if (nIn & nMask)
nOut |= 1;
}
return nOut;
}
namespace {
// Computes and stores twiddle factors for computing DFT later.
struct ScTwiddleFactors
{
ScTwiddleFactors(SCSIZE nN, bool bInverse) :
mfWReal(nN),
mfWImag(nN),
mnN(nN),
mbInverse(bInverse)
{}
void Compute();
void Conjugate()
{
mbInverse = !mbInverse;
for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
mfWImag[nIdx] = -mfWImag[nIdx];
}
std::vector<double> mfWReal;
std::vector<double> mfWImag;
SCSIZE mnN;
bool mbInverse;
};
}
void ScTwiddleFactors::Compute()
{
mfWReal.resize(mnN);
mfWImag.resize(mnN);
double nW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnN);
if (mnN == 1)
{
mfWReal[0] = 1.0;
mfWImag[0] = 0.0;
}
else if (mnN == 2)
{
mfWReal[0] = 1;
mfWImag[0] = 0;
mfWReal[1] = -1;
mfWImag[1] = 0;
}
else if (mnN == 4)
{
mfWReal[0] = 1;
mfWImag[0] = 0;
mfWReal[1] = 0;
mfWImag[1] = (mbInverse ? 1.0 : -1.0);
mfWReal[2] = -1;
mfWImag[2] = 0;
mfWReal[3] = 0;
mfWImag[3] = (mbInverse ? -1.0 : 1.0);
}
else if ((mnN % 4) == 0)
{
const SCSIZE nQSize = mnN >> 2;
// Compute cos of the start quadrant.
// This is the first quadrant if mbInverse == true, else it is the fourth quadrant.
for (SCSIZE nIdx = 0; nIdx <= nQSize; ++nIdx)
mfWReal[nIdx] = cos(nW*static_cast<double>(nIdx));
if (mbInverse)
{
const SCSIZE nQ1End = nQSize;
// First quadrant
for (SCSIZE nIdx = 0; nIdx <= nQ1End; ++nIdx)
mfWImag[nIdx] = mfWReal[nQ1End-nIdx];
// Second quadrant
const SCSIZE nQ2End = nQ1End << 1;
for (SCSIZE nIdx = nQ1End+1; nIdx <= nQ2End; ++nIdx)
{
mfWReal[nIdx] = -mfWReal[nQ2End - nIdx];
mfWImag[nIdx] = mfWImag[nQ2End - nIdx];
}
// Third quadrant
const SCSIZE nQ3End = nQ2End + nQ1End;
for (SCSIZE nIdx = nQ2End+1; nIdx <= nQ3End; ++nIdx)
{
mfWReal[nIdx] = -mfWReal[nIdx - nQ2End];
mfWImag[nIdx] = -mfWImag[nIdx - nQ2End];
}
// Fourth Quadrant
for (SCSIZE nIdx = nQ3End+1; nIdx < mnN; ++nIdx)
{
mfWReal[nIdx] = mfWReal[mnN - nIdx];
mfWImag[nIdx] = -mfWImag[mnN - nIdx];
}
}
else
{
const SCSIZE nQ4End = nQSize;
const SCSIZE nQ3End = nQSize << 1;
const SCSIZE nQ2End = nQ3End + nQSize;
// Fourth quadrant.
for (SCSIZE nIdx = 0; nIdx <= nQ4End; ++nIdx)
mfWImag[nIdx] = -mfWReal[nQ4End - nIdx];
// Third quadrant.
for (SCSIZE nIdx = nQ4End+1; nIdx <= nQ3End; ++nIdx)
{
mfWReal[nIdx] = -mfWReal[nQ3End - nIdx];
mfWImag[nIdx] = mfWImag[nQ3End - nIdx];
}
// Second quadrant.
for (SCSIZE nIdx = nQ3End+1; nIdx <= nQ2End; ++nIdx)
{
mfWReal[nIdx] = -mfWReal[nIdx - nQ3End];
mfWImag[nIdx] = -mfWImag[nIdx - nQ3End];
}
// First quadrant.
for (SCSIZE nIdx = nQ2End+1; nIdx < mnN; ++nIdx)
{
mfWReal[nIdx] = mfWReal[mnN - nIdx];
mfWImag[nIdx] = -mfWImag[mnN - nIdx];
}
}
}
else
{
for (SCSIZE nIdx = 0; nIdx < mnN; ++nIdx)
{
double fAngle = nW*static_cast<double>(nIdx);
mfWReal[nIdx] = cos(fAngle);
mfWImag[nIdx] = sin(fAngle);
}
}
}
namespace {
// A radix-2 decimation in time FFT algorithm for complex valued input.
class ScComplexFFT2
{
public:
// rfArray.size() would always be even and a power of 2. (asserted in prepare())
// rfArray's first half contains the real parts and the later half contains the imaginary parts.
ScComplexFFT2(std::vector<double>& raArray, bool bInverse, bool bPolar, double fMinMag,
ScTwiddleFactors& rTF, bool bSubSampleTFs = false, bool bDisableNormalize = false) :
mrArray(raArray),
mfWReal(rTF.mfWReal),
mfWImag(rTF.mfWImag),
mnPoints(raArray.size()/2),
mnStages(0),
mfMinMag(fMinMag),
mbInverse(bInverse),
mbPolar(bPolar),
mbDisableNormalize(bDisableNormalize),
mbSubSampleTFs(bSubSampleTFs)
{}
void Compute();
private:
void prepare();
double getReal(SCSIZE nIdx)
{
return mrArray[nIdx];
}
void setReal(double fVal, SCSIZE nIdx)
{
mrArray[nIdx] = fVal;
}
double getImag(SCSIZE nIdx)
{
return mrArray[mnPoints + nIdx];
}
void setImag(double fVal, SCSIZE nIdx)
{
mrArray[mnPoints + nIdx] = fVal;
}
SCSIZE getTFactorIndex(SCSIZE nPtIndex, SCSIZE nTfIdxScaleBits)
{
return ( ( nPtIndex << nTfIdxScaleBits ) & ( mnPoints - 1 ) ); // (x & (N-1)) is same as (x % N) but faster.
}
void computeFly(SCSIZE nTopIdx, SCSIZE nBottomIdx, SCSIZE nWIdx1, SCSIZE nWIdx2)
{
if (mbSubSampleTFs)
{
nWIdx1 <<= 1;
nWIdx2 <<= 1;
}
const double x1r = getReal(nTopIdx);
const double x2r = getReal(nBottomIdx);
const double& w1r = mfWReal[nWIdx1];
const double& w1i = mfWImag[nWIdx1];
const double& w2r = mfWReal[nWIdx2];
const double& w2i = mfWImag[nWIdx2];
const double x1i = getImag(nTopIdx);
const double x2i = getImag(nBottomIdx);
setReal(x1r + x2r*w1r - x2i*w1i, nTopIdx);
setImag(x1i + x2i*w1r + x2r*w1i, nTopIdx);
setReal(x1r + x2r*w2r - x2i*w2i, nBottomIdx);
setImag(x1i + x2i*w2r + x2r*w2i, nBottomIdx);
}
std::vector<double>& mrArray;
std::vector<double>& mfWReal;
std::vector<double>& mfWImag;
SCSIZE mnPoints;
SCSIZE mnStages;
double mfMinMag;
bool mbInverse:1;
bool mbPolar:1;
bool mbDisableNormalize:1;
bool mbSubSampleTFs:1;
};
}
void ScComplexFFT2::prepare()
{
SCSIZE nPoints = mnPoints;
lcl_roundUpNearestPow2(nPoints, mnStages);
assert(nPoints == mnPoints);
// Reorder array by bit-reversed indices.
for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
{
SCSIZE nRevIdx = lcl_bitReverse(nIdx, mnPoints);
if (nIdx < nRevIdx)
{
double fTmp = getReal(nIdx);
setReal(getReal(nRevIdx), nIdx);
setReal(fTmp, nRevIdx);
fTmp = getImag(nIdx);
setImag(getImag(nRevIdx), nIdx);
setImag(fTmp, nRevIdx);
}
}
}
static void lcl_normalize(std::vector<double>& rCmplxArray, bool bScaleOnlyReal)
{
const SCSIZE nPoints = rCmplxArray.size()/2;
const double fScale = 1.0/static_cast<double>(nPoints);
// Scale the real part
for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
rCmplxArray[nIdx] *= fScale;
if (!bScaleOnlyReal)
{
const SCSIZE nLen = nPoints*2;
for (SCSIZE nIdx = nPoints; nIdx < nLen; ++nIdx)
rCmplxArray[nIdx] *= fScale;
}
}
static void lcl_convertToPolar(std::vector<double>& rCmplxArray, double fMinMag)
{
const SCSIZE nPoints = rCmplxArray.size()/2;
double fMag, fPhase, fR, fI;
for (SCSIZE nIdx = 0; nIdx < nPoints; ++nIdx)
{
fR = rCmplxArray[nIdx];
fI = rCmplxArray[nPoints+nIdx];
fMag = std::hypot(fR, fI);
if (fMag < fMinMag)
{
fMag = 0.0;
fPhase = 0.0;
}
else
{
fPhase = atan2(fI, fR);
}
rCmplxArray[nIdx] = fMag;
rCmplxArray[nPoints+nIdx] = fPhase;
}
}
void ScComplexFFT2::Compute()
{
prepare();
const SCSIZE nFliesInStage = mnPoints/2;
for (SCSIZE nStage = 0; nStage < mnStages; ++nStage)
{
const SCSIZE nTFIdxScaleBits = mnStages - nStage - 1; // Twiddle factor index's scale factor in bits.
const SCSIZE nFliesInGroup = SCSIZE(1) << nStage;
const SCSIZE nGroups = nFliesInStage/nFliesInGroup;
const SCSIZE nFlyWidth = nFliesInGroup;
for (SCSIZE nGroup = 0, nFlyTopIdx = 0; nGroup < nGroups; ++nGroup)
{
for (SCSIZE nFly = 0; nFly < nFliesInGroup; ++nFly, ++nFlyTopIdx)
{
SCSIZE nFlyBottomIdx = nFlyTopIdx + nFlyWidth;
SCSIZE nWIdx1 = getTFactorIndex(nFlyTopIdx, nTFIdxScaleBits);
SCSIZE nWIdx2 = getTFactorIndex(nFlyBottomIdx, nTFIdxScaleBits);
computeFly(nFlyTopIdx, nFlyBottomIdx, nWIdx1, nWIdx2);
}
nFlyTopIdx += nFlyWidth;
}
}
if (mbPolar)
lcl_convertToPolar(mrArray, mfMinMag);
// Normalize after converting to polar, so we have a chance to
// save O(mnPoints) flops.
if (mbInverse && !mbDisableNormalize)
lcl_normalize(mrArray, mbPolar);
}
namespace {
// Bluestein's algorithm or chirp z-transform algorithm that can be used to
// compute DFT of a complex valued input of any length N in O(N lgN) time.
class ScComplexBluesteinFFT
{
public:
ScComplexBluesteinFFT(std::vector<double>& rArray, bool bReal, bool bInverse,
bool bPolar, double fMinMag, bool bDisableNormalize = false) :
mrArray(rArray),
mnPoints(rArray.size()/2), // rArray should have space for imaginary parts even if real input.
mfMinMag(fMinMag),
mbReal(bReal),
mbInverse(bInverse),
mbPolar(bPolar),
mbDisableNormalize(bDisableNormalize)
{}
void Compute();
private:
std::vector<double>& mrArray;
const SCSIZE mnPoints;
double mfMinMag;
bool mbReal:1;
bool mbInverse:1;
bool mbPolar:1;
bool mbDisableNormalize:1;
};
}
void ScComplexBluesteinFFT::Compute()
{
std::vector<double> aRealScalars(mnPoints);
std::vector<double> aImagScalars(mnPoints);
double fW = (mbInverse ? 2 : -2)*M_PI/static_cast<double>(mnPoints);
for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
{
double fAngle = 0.5*fW*static_cast<double>(nIdx*nIdx);
aRealScalars[nIdx] = cos(fAngle);
aImagScalars[nIdx] = sin(fAngle);
}
SCSIZE nMinSize = mnPoints*2 - 1;
SCSIZE nExtendedLength = nMinSize, nTmp = 0;
lcl_roundUpNearestPow2(nExtendedLength, nTmp);
std::vector<double> aASignal(nExtendedLength*2); // complex valued
std::vector<double> aBSignal(nExtendedLength*2); // complex valued
double fReal, fImag;
for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
{
// Real part of A signal.
aASignal[nIdx] = mrArray[nIdx]*aRealScalars[nIdx] + (mbReal ? 0.0 : -mrArray[mnPoints+nIdx]*aImagScalars[nIdx]);
// Imaginary part of A signal.
aASignal[nExtendedLength + nIdx] = mrArray[nIdx]*aImagScalars[nIdx] + (mbReal ? 0.0 : mrArray[mnPoints+nIdx]*aRealScalars[nIdx]);
// Real part of B signal.
aBSignal[nIdx] = fReal = aRealScalars[nIdx];
// Imaginary part of B signal.
aBSignal[nExtendedLength + nIdx] = fImag = -aImagScalars[nIdx]; // negative sign because B signal is the conjugation of the scalars.
if (nIdx)
{
// B signal needs a mirror of its part in 0 < n < mnPoints at the tail end.
aBSignal[nExtendedLength - nIdx] = fReal;
aBSignal[(nExtendedLength<<1) - nIdx] = fImag;
}
}
{
ScTwiddleFactors aTF(nExtendedLength, false /*not inverse*/);
aTF.Compute();
// Do complex-FFT2 of both A and B signal.
ScComplexFFT2 aFFT2A(aASignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
aTF, false /*no subsample*/, true /* disable normalize */);
aFFT2A.Compute();
ScComplexFFT2 aFFT2B(aBSignal, false /*not inverse*/, false /*no polar*/, 0.0 /* no clipping */,
aTF, false /*no subsample*/, true /* disable normalize */);
aFFT2B.Compute();
double fAR, fAI, fBR, fBI;
for (SCSIZE nIdx = 0; nIdx < nExtendedLength; ++nIdx)
{
fAR = aASignal[nIdx];
fAI = aASignal[nExtendedLength + nIdx];
fBR = aBSignal[nIdx];
fBI = aBSignal[nExtendedLength + nIdx];
// Do point-wise product inplace in A signal.
aASignal[nIdx] = fAR*fBR - fAI*fBI;
aASignal[nExtendedLength + nIdx] = fAR*fBI + fAI*fBR;
}
// Do complex-inverse-FFT2 of aASignal.
aTF.Conjugate();
ScComplexFFT2 aFFT2AI(aASignal, true /*inverse*/, false /*no polar*/, 0.0 /* no clipping */, aTF); // Need normalization here.
aFFT2AI.Compute();
}
// Point-wise multiply with scalars.
for (SCSIZE nIdx = 0; nIdx < mnPoints; ++nIdx)
{
fReal = aASignal[nIdx];
fImag = aASignal[nExtendedLength + nIdx];
mrArray[nIdx] = fReal*aRealScalars[nIdx] - fImag*aImagScalars[nIdx]; // no conjugation needed here.
mrArray[mnPoints + nIdx] = fReal*aImagScalars[nIdx] + fImag*aRealScalars[nIdx];
}
// Normalize/Polar operations
if (mbPolar)
lcl_convertToPolar(mrArray, mfMinMag);
// Normalize after converting to polar, so we have a chance to
// save O(mnPoints) flops.
if (mbInverse && !mbDisableNormalize)
lcl_normalize(mrArray, mbPolar);
}
namespace {
// Computes DFT of an even length(N) real-valued input by using a
// ScComplexFFT2 if N == 2^k for some k or else by using a ScComplexBluesteinFFT
// with a complex valued input of length = N/2.
class ScRealFFT
{
public:
ScRealFFT(std::vector<double>& rInArray, std::vector<double>& rOutArray, bool bInverse,
bool bPolar, double fMinMag) :
mrInArray(rInArray),
mrOutArray(rOutArray),
mfMinMag(fMinMag),
mbInverse(bInverse),
mbPolar(bPolar)
{}
void Compute();
private:
std::vector<double>& mrInArray;
std::vector<double>& mrOutArray;
double mfMinMag;
bool mbInverse:1;
bool mbPolar:1;
};
}
void ScRealFFT::Compute()
{
// input length has to be even to do this optimization.
assert(mrInArray.size() % 2 == 0);
assert(mrInArray.size()*2 == mrOutArray.size());
// nN is the number of points in the complex-fft input
// which will be half of the number of points in real array.
const SCSIZE nN = mrInArray.size()/2;
if (nN == 0)
{
mrOutArray[0] = mrInArray[0];
mrOutArray[1] = 0.0;
return;
}
// work array should be the same length as mrInArray
std::vector<double> aWorkArray(nN*2);
for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
{
SCSIZE nDoubleIdx = 2*nIdx;
// Use even elements as real part
aWorkArray[nIdx] = mrInArray[nDoubleIdx];
// and odd elements as imaginary part of the contrived complex sequence.
aWorkArray[nN+nIdx] = mrInArray[nDoubleIdx+1];
}
ScTwiddleFactors aTFs(nN*2, mbInverse);
aTFs.Compute();
SCSIZE nNextPow2 = nN, nTmp = 0;
lcl_roundUpNearestPow2(nNextPow2, nTmp);
if (nNextPow2 == nN)
{
ScComplexFFT2 aFFT2(aWorkArray, mbInverse, false /*disable polar*/, 0.0 /* no clipping */,
aTFs, true /*subsample tf*/, true /*disable normalize*/);
aFFT2.Compute();
}
else
{
ScComplexBluesteinFFT aFFT(aWorkArray, false /*complex input*/, mbInverse, false /*disable polar*/,
0.0 /* no clipping */, true /*disable normalize*/);
aFFT.Compute();
}
// Post process aWorkArray to populate mrOutArray
const SCSIZE nTwoN = 2*nN, nThreeN = 3*nN;
double fY1R, fY2R, fY1I, fY2I, fResR, fResI, fWR, fWI;
for (SCSIZE nIdx = 0; nIdx < nN; ++nIdx)
{
const SCSIZE nIdxRev = nIdx ? (nN - nIdx) : 0;
fY1R = aWorkArray[nIdx];
fY2R = aWorkArray[nIdxRev];
fY1I = aWorkArray[nN + nIdx];
fY2I = aWorkArray[nN + nIdxRev];
fWR = aTFs.mfWReal[nIdx];
fWI = aTFs.mfWImag[nIdx];
// mrOutArray has length = 4*nN
// Real part of the final output (only half of the symmetry around Nyquist frequency)
// Fills the first quarter.
mrOutArray[nIdx] = fResR = 0.5*(
fY1R + fY2R +
fWR * (fY1I + fY2I) +
fWI * (fY1R - fY2R) );
// Imaginary part of the final output (only half of the symmetry around Nyquist frequency)
// Fills the third quarter.
mrOutArray[nTwoN + nIdx] = fResI = 0.5*(
fY1I - fY2I +
fWI * (fY1I + fY2I) -
fWR * (fY1R - fY2R) );
// Fill the missing 2 quarters using symmetry argument.
if (nIdx)
{
// Fills the 2nd quarter.
mrOutArray[nN + nIdxRev] = fResR;
// Fills the 4th quarter.
mrOutArray[nThreeN + nIdxRev] = -fResI;
}
else
{
mrOutArray[nN] = fY1R - fY1I;
mrOutArray[nThreeN] = 0.0;
}
}
// Normalize/Polar operations
if (mbPolar)
lcl_convertToPolar(mrOutArray, mfMinMag);
// Normalize after converting to polar, so we have a chance to
// save O(mnPoints) flops.
if (mbInverse)
lcl_normalize(mrOutArray, mbPolar);
}
using ScMatrixGenerator = ScMatrixRef(SCSIZE, SCSIZE, std::vector<double>&);
namespace {
// Generic FFT class that decides which FFT implementation to use.
class ScFFT
{
public:
ScFFT(ScMatrixRef& pMat, bool bReal, bool bInverse, bool bPolar, double fMinMag) :
mpInputMat(pMat),
mfMinMag(fMinMag),
mbReal(bReal),
mbInverse(bInverse),
mbPolar(bPolar)
{}
ScMatrixRef Compute(const std::function<ScMatrixGenerator>& rMatGenFunc);
private:
ScMatrixRef& mpInputMat;
double mfMinMag;
bool mbReal:1;
bool mbInverse:1;
bool mbPolar:1;
};
}
ScMatrixRef ScFFT::Compute(const std::function<ScMatrixGenerator>& rMatGenFunc)
{
std::vector<double> aArray;
mpInputMat->GetDoubleArray(aArray);
SCSIZE nPoints = mbReal ? aArray.size() : (aArray.size()/2);
if (nPoints == 1)
{
std::vector<double> aOutArray(2);
aOutArray[0] = aArray[0];
aOutArray[1] = mbReal ? 0.0 : aArray[1];
if (mbPolar)
lcl_convertToPolar(aOutArray, mfMinMag);
return rMatGenFunc(2, 1, aOutArray);
}
if (mbReal && (nPoints % 2) == 0)
{
std::vector<double> aOutArray(nPoints*2);
ScRealFFT aFFT(aArray, aOutArray, mbInverse, mbPolar, mfMinMag);
aFFT.Compute();
return rMatGenFunc(2, nPoints, aOutArray);
}
SCSIZE nNextPow2 = nPoints, nTmp = 0;
lcl_roundUpNearestPow2(nNextPow2, nTmp);
if (nNextPow2 == nPoints && !mbReal)
{
ScTwiddleFactors aTF(nPoints, mbInverse);
aTF.Compute();
ScComplexFFT2 aFFT2(aArray, mbInverse, mbPolar, mfMinMag, aTF);
aFFT2.Compute();
return rMatGenFunc(2, nPoints, aArray);
}
if (mbReal)
aArray.resize(nPoints*2, 0.0);
ScComplexBluesteinFFT aFFT(aArray, mbReal, mbInverse, mbPolar, mfMinMag);
aFFT.Compute();
return rMatGenFunc(2, nPoints, aArray);
}
void ScInterpreter::ScFourier()
{
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 2, 5 ) )
return;
bool bInverse = false;
bool bPolar = false;
double fMinMag = 0.0;
if (nParamCount == 5)
{
if (IsMissing())
Pop();
else
fMinMag = GetDouble();
}
if (nParamCount >= 4)
{
if (IsMissing())
Pop();
else
bPolar = GetBool();
}
if (nParamCount >= 3)
{
if (IsMissing())
Pop();
else
bInverse = GetBool();
}
bool bGroupedByColumn = GetBool();
ScMatrixRef pInputMat = GetMatrix();
if (!pInputMat)
{
PushIllegalParameter();
return;
}
SCSIZE nC, nR;
pInputMat->GetDimensions(nC, nR);
if ((bGroupedByColumn && nC > 2) || (!bGroupedByColumn && nR > 2))
{
// There can be no more than 2 columns (real, imaginary) if data grouped by columns.
// and no more than 2 rows if data is grouped by rows.
PushIllegalArgument();
return;
}
if (!pInputMat->IsNumeric())
{
PushNoValue();
return;
}
bool bRealInput = true;
if (!bGroupedByColumn)
{
pInputMat->MatTrans(*pInputMat);
bRealInput = (nR == 1);
}
else
{
bRealInput = (nC == 1);
}
ScFFT aFFT(pInputMat, bRealInput, bInverse, bPolar, fMinMag);
std::function<ScMatrixGenerator> aFunc = [this](SCSIZE nCol, SCSIZE nRow, std::vector<double>& rVec) -> ScMatrixRef
{
return this->GetNewMat(nCol, nRow, rVec);
};
ScMatrixRef pOut = aFFT.Compute(aFunc);
PushMatrix(pOut);
}
/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
View git blame Copy permalink