/* -*- 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 .
*/
// The idea how this group of gamma functions is calculated, is // based on the Cephes library // online http://www.moshier.net/#Cephes [called 2008-02]
/** You must ensure fA>0.0 && fX>0.0 valid results only if fX > fA+1.0
uses continued fraction with odd items */ double ScInterpreter::GetGammaContFraction( double fA, double fX )
{
/** You must ensure fA>0.0 && fX>0.0 valid results only if fX <= fA+1.0
uses power series */ double ScInterpreter::GetGammaSeries( double fA, double fX )
{ double fDenomfactor = fA; double fSummand = 1.0/fA; double fSum = fSummand; int nCount=1; do
{
fDenomfactor = fDenomfactor + 1.0;
fSummand = fSummand * fX/fDenomfactor;
fSum = fSum + fSummand;
nCount = nCount+1;
} while ( fSummand/fSum > fHalfMachEps && nCount<=10000); // large amount of iterations will be carried out for huge fAlpha, even // if fX <= fAlpha+1.0 if (nCount>10000)
{
SetError(FormulaError::NoConvergence);
} return fSum;
}
/** You must ensure fA>0.0 && fX>0.0) */ double ScInterpreter::GetLowRegIGamma( double fA, double fX )
{ double fLnFactor = fA * log(fX) - fX - GetLogGamma(fA); double fFactor = exp(fLnFactor); // Do we need more accuracy than exp(ln()) has? if (fX>fA+1.0) // includes fX>1.0; 1-GetUpRegIGamma, continued fraction return 1.0 - fFactor * GetGammaContFraction(fA,fX); else// fX<=1.0 || fX<=fA+1.0, series return fFactor * GetGammaSeries(fA,fX);
}
/** You must ensure fA>0.0 && fX>0.0) */ double ScInterpreter::GetUpRegIGamma( double fA, double fX )
{
double fLnFactor= fA*log(fX)-fX-GetLogGamma(fA); double fFactor = exp(fLnFactor); //Do I need more accuracy than exp(ln()) has?; if (fX>fA+1.0) // includes fX>1.0 return fFactor * GetGammaContFraction(fA,fX); else//fX<=1 || fX<=fA+1, 1-GetLowRegIGamma, series return 1.0 -fFactor * GetGammaSeries(fA,fX);
}
/** Gamma distribution, probability density function. fLambda is "scale" parameter
You must ensure fAlpha>0.0 and fLambda>0.0 */ double ScInterpreter::GetGammaDistPDF( double fX, double fAlpha, double fLambda )
{ if (fX < 0.0) return 0.0; // see ODFF elseif (fX == 0) // in this case 0^0 isn't zero
{ if (fAlpha < 1.0)
{
SetError(FormulaError::DivisionByZero); // should be #DIV/0 return HUGE_VAL;
} elseif (fAlpha == 1)
{ return (1.0 / fLambda);
} else
{ return 0.0;
}
} else
{ double fXr = fX / fLambda; // use exp(ln()) only for large arguments because of less accuracy if (fXr > 1.0)
{ constdouble fLogDblMax = log( ::std::numeric_limits<double>::max()); if (log(fXr) * (fAlpha-1.0) < fLogDblMax && fAlpha < fMaxGammaArgument)
{ return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha);
} else
{ return exp( (fAlpha-1.0) * log(fXr) - fXr - log(fLambda) - GetLogGamma(fAlpha));
}
} else// fXr near to zero
{ if (fAlpha<fMaxGammaArgument)
{ return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha);
} else
{ return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / exp( GetLogGamma(fAlpha));
}
}
}
}
/** Gamma distribution, cumulative distribution function. fLambda is "scale" parameter
You must ensure fAlpha>0.0 and fLambda>0.0 */ double ScInterpreter::GetGammaDist( double fX, double fAlpha, double fLambda )
{ if (fX <= 0.0) return 0.0; else return GetLowRegIGamma( fAlpha, fX / fLambda);
}
namespace {
class NumericCellAccumulator
{
KahanSum maSum;
FormulaError mnError;
// Get the number format of the last iterated cell.
nFuncFmtIndex = aAction.getNumberFormat();
} else
{
FuncCount aAction(mrContext);
aSet.executeColumnAction(mrDoc, aAction);
nCount += aAction.getCount();
// Get the number format of the last iterated cell.
nFuncFmtIndex = aAction.getNumberFormat();
}
nFuncFmtType = mrContext.NFGetType(nFuncFmtIndex);
} else
{
ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags, bTextAsZero );
FormulaError nErr = FormulaError::NONE; if (aValIter.GetFirst(fVal, nErr))
{ // placed the loop on the inside for performance reasons:
aValIter.GetCurNumFmtInfo( nFuncFmtType, nFuncFmtIndex ); switch( eFunc )
{ case ifAVERAGE: case ifSUM: if ( mnSubTotalFlags & SubtotalFlags::IgnoreErrVal )
{ do
{ if ( nErr == FormulaError::NONE )
{
SetError(nErr);
fRes += fVal;
nCount++;
}
} while (aValIter.GetNext(fVal, nErr));
} else
{ do
{
SetError(nErr);
fRes += fVal;
nCount++;
} while (aValIter.GetNext(fVal, nErr));
} break; case ifSUMSQ: if ( mnSubTotalFlags & SubtotalFlags::IgnoreErrVal )
{ do
{ if ( nErr == FormulaError::NONE )
{
SetError(nErr);
fRes += fVal * fVal;
nCount++;
}
} while (aValIter.GetNext(fVal, nErr));
} else
{ do
{
SetError(nErr);
fRes += fVal * fVal;
nCount++;
} while (aValIter.GetNext(fVal, nErr));
} break; case ifPRODUCT: do
{ if ( !( nErr != FormulaError::NONE && ( mnSubTotalFlags & SubtotalFlags::IgnoreErrVal ) ) )
{
SetError(nErr);
fRes *= fVal;
nCount++;
}
} while (aValIter.GetNext(fVal, nErr)); break; case ifCOUNT: do
{ if ( nErr == FormulaError::NONE )
nCount++;
} while (aValIter.GetNext(fVal, nErr)); break; default: ; // nothing
}
SetError( nErr );
}
} if (nRefArrayPos != std::numeric_limits<size_t>::max())
{ // Update vector element with current value. if (xResCount)
xResCount->PutDouble( xResCount->GetDouble(0,nRefArrayPos) + nCount, 0,nRefArrayPos); double fVecRes = xResMat->GetDouble(0,nRefArrayPos); if (eFunc == ifPRODUCT)
fVecRes *= fRes.get(); else
fVecRes += fRes.get();
xResMat->PutDouble( fVecRes, 0,nRefArrayPos); // Reset.
fRes = ResInitVal;
nCount = 0;
nRefArrayPos = std::numeric_limits<size_t>::max();
}
} break; case svExternalDoubleRef:
{
ScMatrixRef pMat;
PopExternalDoubleRef(pMat); if ( nGlobalError != FormulaError::NONE && !( mnSubTotalFlags & SubtotalFlags::IgnoreErrVal ) ) break;
// A boolean return type makes no sense on sums et al. // Counts are always numbers. if( nFuncFmtType == SvNumFormatType::LOGICAL || eFunc == ifCOUNT || eFunc == ifCOUNT2 )
nFuncFmtType = SvNumFormatType::NUMBER;
if (xResMat)
{ // Include value of last non-references-array type and calculate final result. for (SCSIZE i=0; i < nMatRows; ++i)
{
sal_uLong nVecCount = (xResCount ? nCount + xResCount->GetDouble(0,i) : nCount); double fVecRes = xResMat->GetDouble(0,i); if (eFunc == ifPRODUCT)
fVecRes *= fRes.get(); else
fVecRes += fRes.get();
fVecRes = lcl_IterResult( eFunc, fVecRes, nVecCount);
xResMat->PutDouble( fVecRes, 0,i);
}
PushMatrix( xResMat);
} else
{
PushDouble( lcl_IterResult( eFunc, fRes.get(), nCount));
}
}
/** * The purpose of RAWSUBTRACT() is exactly to not apply any error correction, approximation etc. * But use the "raw" IEEE 754 double subtraction. * So no Kahan summation
*/ void ScInterpreter::ScRawSubtract()
{ short nParamCount = GetByte(); if (!MustHaveParamCountMin( nParamCount, 2)) return;
// Reverse stack to process arguments from left to right.
ReverseStack( nParamCount); // Obtain the minuend. double fRes = GetDouble();
while (nGlobalError == FormulaError::NONE && --nParamCount > 0)
{ // Simple single values without matrix support.
fRes -= GetDouble();
} while (nParamCount-- > 0)
PopError();
