summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorRĂ¼diger Timm <rt@openoffice.org>2008-07-08 06:21:03 +0000
committerRĂ¼diger Timm <rt@openoffice.org>2008-07-08 06:21:03 +0000
commit07e6b192bc9e13ecf29bffda40ef371fc4a944a0 (patch)
treed28187cacd27d9aaff4c5db32acc4696d2ce512a
parent97f95c196906f276a5f06fd6a4bab8411a257f3a (diff)
INTEGRATION: CWS odff04 (1.23.2); FILE MERGED
2008/06/23 10:46:15 er 1.23.2.2: MSVC dumbness, cannot initialize static const non-integral inside class 2008/06/18 15:48:57 er 1.23.2.1: #i90703# accuracy of GAMMADIST, CHIDIST, GAMMAINV and CHIINV; patch from <regina>
-rw-r--r--sc/source/core/tool/interpr3.cxx433
1 files changed, 293 insertions, 140 deletions
diff --git a/sc/source/core/tool/interpr3.cxx b/sc/source/core/tool/interpr3.cxx
index eaaa7f0f1bba..624d3efdae88 100644
--- a/sc/source/core/tool/interpr3.cxx
+++ b/sc/source/core/tool/interpr3.cxx
@@ -7,7 +7,7 @@
* OpenOffice.org - a multi-platform office productivity suite
*
* $RCSfile: interpr3.cxx,v $
- * $Revision: 1.23 $
+ * $Revision: 1.24 $
*
* This file is part of OpenOffice.org.
*
@@ -61,6 +61,8 @@ using ::std::vector;
// PI jetzt als F_PI aus solar.h
//#define PI 3.1415926535897932
+const double ScInterpreter::fMaxGammaArgument = 171.624376956302; // found experimental
+
//-----------------------------------------------------------------------------
class ScDistFunc
@@ -71,94 +73,111 @@ public:
// iteration for inverse distributions
-//template< class T > double lcl_IterateInverse( const T& rFunction, double x0, double x1, BOOL& rConvError )
-double lcl_IterateInverse( const ScDistFunc& rFunction, double x0, double x1, BOOL& rConvError )
+//template< class T > double lcl_IterateInverse( const T& rFunction, double x0, double x1, bool& rConvError )
+
+/** u*w<0.0 fails for values near zero */
+inline bool lcl_HasChangeOfSign( double u, double w )
+{
+ return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);
+}
+
+double lcl_IterateInverse( const ScDistFunc& rFunction, double fAx, double fBx, bool& rConvError )
{
- rConvError = FALSE;
- double fEps = 1.0E-7;
+ rConvError = false;
+ const double fYEps = 1.0E-307;
+ const double fXEps = ::std::numeric_limits<double>::epsilon();
- DBG_ASSERT(x0<x1, "IterateInverse: wrong interval");
+ DBG_ASSERT(fAx<fBx, "IterateInverse: wrong interval");
// find enclosing interval
- double f0 = rFunction.GetValue(x0);
- double f1 = rFunction.GetValue(x1);
- double xs;
- USHORT i;
- for (i = 0; i < 1000 && f0*f1 > 0.0; i++)
+ double fAy = rFunction.GetValue(fAx);
+ double fBy = rFunction.GetValue(fBx);
+ double fTemp;
+ unsigned short nCount;
+ for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy); nCount++)
{
- if (fabs(f0) <= fabs(f1))
+ if (fabs(fAy) <= fabs(fBy))
{
- xs = x0;
- x0 += 2.0 * (x0 - x1);
- if (x0 < 0.0)
- x0 = 0.0;
- x1 = xs;
- f1 = f0;
- f0 = rFunction.GetValue(x0);
+ fTemp = fAx;
+ fAx += 2.0 * (fAx - fBx);
+ if (fAx < 0.0)
+ fAx = 0.0;
+ fBx = fTemp;
+ fBy = fAy;
+ fAy = rFunction.GetValue(fAx);
}
else
{
- xs = x1;
- x1 += 2.0 * (x1 - x0);
- x0 = xs;
- f0 = f1;
- f1 = rFunction.GetValue(x1);
+ fTemp = fBx;
+ fBx += 2.0 * (fBx - fAx);
+ fAx = fTemp;
+ fAy = fBy;
+ fBy = rFunction.GetValue(fBx);
}
}
- if (f0 == 0.0)
- return x0;
- if (f1 == 0.0)
- return x1;
-
- // simple iteration
-
- double x00 = x0;
- double x11 = x1;
- double fs = 0.0;
- for (i = 0; i < 100; i++)
+ if (fAy == 0.0)
+ return fAx;
+ if (fBy == 0.0)
+ return fBx;
+ if (!lcl_HasChangeOfSign( fAy, fBy))
{
- xs = 0.5*(x0+x1);
- if (fabs(f1-f0) >= fEps)
+ rConvError = true;
+ return 0.0;
+ }
+ // inverse quadric interpolation with additional brackets
+ // set three points
+ double fPx = fAx;
+ double fPy = fAy;
+ double fQx = fBx;
+ double fQy = fBy;
+ double fRx = fAx;
+ double fRy = fAy;
+ double fSx = 0.5 * (fAx + fBx); // potential next point
+ bool bHasToInterpolate = true;
+ nCount = 0;
+ while ( nCount < 500 && fabs(fRy) > fYEps &&
+ (fBx-fAx) > ::std::max( fabs(fAx), fabs(fBx)) * fXEps )
+ {
+ if (bHasToInterpolate)
{
- fs = rFunction.GetValue(xs);
- if (f0*fs <= 0.0)
+ if (fPy!=fQy && fQy!=fRy && fRy!=fPy)
{
- x1 = xs;
- f1 = fs;
+ 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
- {
- x0 = xs;
- f0 = fs;
- }
+ bHasToInterpolate = false;
+ }
+ if(!bHasToInterpolate)
+ {
+ fSx = 0.5 * (fAx + fBx);
+ // reset points
+ fPx = fAx; fPy = fAy;
+ fQx = fBx; fQy = fBy;
+ bHasToInterpolate = true;
+ }
+ // shift points for next interpolation
+ fPx = fQx; fQx = fRx; fRx = fSx;
+ fPy = fQy; fQy = fRy; fRy = rFunction.GetValue(fSx);
+ // update brackets
+ if (lcl_HasChangeOfSign( fAy, fRy))
+ {
+ fBx = fRx; fBy = fRy;
}
else
{
- // add one step of regula falsi to improve precision
-
- if ( x0 != x1 )
- {
- double regxs = (f1-f0)/(x1-x0);
- if ( regxs != 0.0)
- {
- double regx = x1 - f1/regxs;
- if (regx >= x00 && regx <= x11)
- {
- double regfs = rFunction.GetValue(regx);
- if ( fabs(regfs) < fabs(fs) )
- xs = regx;
- }
- }
- }
-
- return xs;
+ fAx = fRx; fAy = fRy;
}
+ // if last interration brought to small advance, then do bisection next
+ // time, for safety
+ bHasToInterpolate = bHasToInterpolate && (fabs(fRy) * 2.0 <= fabs(fQy));
+ ++nCount;
}
-
- rConvError = TRUE;
- return 0.0;
+ return fRx;
}
//-----------------------------------------------------------------------------
@@ -475,50 +494,160 @@ double ScInterpreter::BinomKoeff(double n, double k)
return nVal;
}
-double ScInterpreter::GammaHelp(double& x, BOOL& bReflect)
+
+// The algorithm is based on lanczos13m53 in lanczos.hpp
+// in math library from http://www.boost.org
+/** you must ensure fZ>0
+ Uses a variant of the Lanczos sum with a rational function. */
+double lcl_getLanczosSum(double fZ)
{
- double c[6] = {76.18009173, -86.50532033, 24.01409822,
- -1.231739516, 0.120858003E-2, -0.536382E-5};
- if (x >= 1.0)
- {
- bReflect = FALSE;
- x -= 1.0;
+ const double fNum[13] ={
+ 23531376880.41075968857200767445163675473,
+ 42919803642.64909876895789904700198885093,
+ 35711959237.35566804944018545154716670596,
+ 17921034426.03720969991975575445893111267,
+ 6039542586.35202800506429164430729792107,
+ 1439720407.311721673663223072794912393972,
+ 248874557.8620541565114603864132294232163,
+ 31426415.58540019438061423162831820536287,
+ 2876370.628935372441225409051620849613599,
+ 186056.2653952234950402949897160456992822,
+ 8071.672002365816210638002902272250613822,
+ 210.8242777515793458725097339207133627117,
+ 2.506628274631000270164908177133837338626
+ };
+ const double fDenom[13] = {
+ 0,
+ 39916800,
+ 120543840,
+ 150917976,
+ 105258076,
+ 45995730,
+ 13339535,
+ 2637558,
+ 357423,
+ 32670,
+ 1925,
+ 66,
+ 1
+ };
+ // Horner scheme
+ double fSumNum;
+ double fSumDenom;
+ int nI;
+ double fZInv;
+ 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
{
- bReflect = TRUE;
- x = 1.0 - x;
- }
- double s, anum;
- s = 1.0;
- anum = x;
- for (USHORT i = 0; i < 6; i++)
- {
- anum += 1.0;
- s += c[i]/anum;
+ 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];
+ }
}
- s *= 2.506628275; // sqrt(2*PI)
- return s;
+ return fSumNum/fSumDenom;
}
-double ScInterpreter::GetGamma(double x)
+// The algorithm is based on tgamma in gamma.hpp
+// in math library from http://www.boost.org
+/** You must ensure fZ>0; fZ>171.624376956302 will overflow. */
+double lcl_GetGammaHelper(double fZ)
{
- BOOL bReflect;
- double G = GammaHelp(x, bReflect);
- G = pow(x+5.5,x+0.5)*G/exp(x+5.5);
- if (bReflect)
- G = F_PI*x/(G*::rtl::math::sin(F_PI*x));
- return G;
+ 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;
}
-double ScInterpreter::GetLogGamma(double x)
+// The algorithm is based on tgamma in gamma.hpp
+// in math library from http://www.boost.org
+/** You must ensure fZ>0 */
+double lcl_GetLogGammaHelper(double fZ)
{
- BOOL bReflect;
- double G = GammaHelp(x, bReflect);
- G = (x+0.5)*log(x+5.5)+log(G)-(x+5.5);
- if (bReflect)
- G = log(F_PI*x)-G-log(::rtl::math::sin(F_PI*x));
- return G;
+ 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(F_PI);
+ const double fLogDblMax = log( ::std::numeric_limits<double>::max());
+
+ if (fZ > fMaxGammaArgument)
+ {
+ SetError(errIllegalFPOperation);
+ return HUGE_VAL;
+ }
+
+ if (fZ >= 1.0)
+ return lcl_GetGammaHelper(fZ);
+
+ if (fZ >= 0.5) // shift to x>=1 using Gamma(x)=Gamma(x+1)/x
+ return lcl_GetGammaHelper(fZ+1) / fZ;
+
+ if (fZ >= -0.5) // shift to x>=1, might overflow
+ {
+ double fLogTest = lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log( fabs(fZ));
+ if (fLogTest >= fLogDblMax)
+ {
+ SetError( errIllegalFPOperation);
+ return HUGE_VAL;
+ }
+ return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ;
+ }
+ // fZ<-0.5
+ // Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)*sin(pi*x) )
+ double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( fabs( ::rtl::math::sin( F_PI*fZ)));
+ if (fLogDivisor - fLogPi >= fLogDblMax) // underflow
+ return 0.0;
+
+ if (fLogDivisor<0.0)
+ if (fLogPi - fLogDivisor > fLogDblMax) // overflow
+ {
+ SetError(errIllegalFPOperation);
+ return HUGE_VAL;
+ }
+
+ return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( F_PI*fZ) < 0.0) ? -1.0 : 1.0);
+}
+
+
+/** You must ensure fZ>0 */
+double ScInterpreter::GetLogGamma(double fZ)
+{
+ if (fZ >= fMaxGammaArgument)
+ return lcl_GetLogGammaHelper(fZ);
+ if (fZ >= 1.0)
+ return log(lcl_GetGammaHelper(fZ));
+ if (fZ >= 0.5)
+ return log( lcl_GetGammaHelper(fZ+1) / fZ);
+ return lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log(fZ);
}
double ScInterpreter::GetBetaDist(double x, double alpha, double beta)
@@ -649,42 +778,58 @@ double ScInterpreter::GetTDist(double T, double fDF)
*/
}
-double ScInterpreter::GetChiDist(double fChi, double fDF)
+// for LEGACY.CHIDIST, returns right tail, fDF=degrees of freedom
+/** You must ensure fDF>0.0 */
+double ScInterpreter::GetChiDist(double fX, double fDF)
{
- return 1.0 - GetGammaDist(fChi/2.0, fDF/2.0, 1.0);
-/*
- double x = 1.0;
- for (double i = fDF; i >= 2.0; i -= 2.0)
- x *= fChi/i;
- x *= exp(-fChi/2.0);
- if (fmod(fDF, 2.0) != 0.0)
- x *= sqrt(2.0*fChi/F_PI);
- double S = 1.0;
- double T = 1.0;
- double G = fDF;
- BOOL bStop = FALSE;
- while (!bStop)
- {
- G += 2.0;
- T *= fChi/G;
- if (T < 1.0E-7)
- bStop = TRUE;
- else
- S += T;
+ 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);
+}
+
+// for ODF 1.2 GAMMA
+void ScInterpreter::ScGamma()
+{
+ double x = GetDouble();
+ double fResult;
+ if (x <= 0.0 && x == ::rtl::math::approxFloor(x))
+ PushIllegalArgument();
+ else
+ {
+ fResult = GetGamma(x);
+ if (nGlobalError)
+ {
+ PushError( nGlobalError);
+ return;
+ }
+ PushDouble(fResult);
}
- return 1.0 - x*S;
-*/
}
+
void ScInterpreter::ScLogGamma()
{
double x = GetDouble();
- if (x > 0.0)
- PushDouble(GetLogGamma(x));
+ if (x > 0.0) // constraint from ODFF
+ PushDouble( GetLogGamma(x));
else
PushIllegalArgument();
}
+
void ScInterpreter::ScBetaDist()
{
BYTE nParamCount = GetByte();
@@ -1176,6 +1321,7 @@ void ScInterpreter::ScFDist()
void ScInterpreter::ScChiDist()
{
+ double fResult;
if ( !MustHaveParamCount( GetByte(), 2 ) )
return;
double fDF = ::rtl::math::approxFloor(GetDouble());
@@ -1185,7 +1331,13 @@ void ScInterpreter::ScChiDist()
PushIllegalArgument();
return;
}
- PushDouble(GetChiDist(fChi, fDF));
+ fResult = GetChiDist( fChi, fDF);
+ if (nGlobalError)
+ {
+ PushError( nGlobalError);
+ return;
+ }
+ PushDouble(fResult);
}
void ScInterpreter::ScWeibull()
@@ -1491,19 +1643,20 @@ void ScInterpreter::ScGammaDist()
{
if ( !MustHaveParamCount( GetByte(), 4 ) )
return;
- double kum = GetDouble(); // 0 oder 1
- double beta = GetDouble();
- double alpha = GetDouble();
- double x = GetDouble(); // x
- if (x < 0.0 || alpha <= 0.0 || beta <= 0.0)
+ double fCumulative = GetDouble(); // parameter type will change to bool
+ bool bCumulative = (fCumulative == 1.0); // see ODFF
+ double fBeta = GetDouble(); // scale
+ double fAlpha = GetDouble(); // shape
+ double fX = GetDouble(); // x
+ if (fAlpha <= 0.0 || fBeta <= 0.0)
PushIllegalArgument();
- else if (kum == 0.0) // Dichte
+ else
{
- double G = GetGamma(alpha);
- PushDouble(pow(x,alpha-1.0)/exp(x/beta)/pow(beta,alpha)/G);
+ if (bCumulative) // distribution
+ PushDouble( GetGammaDist( fX, fAlpha, fBeta));
+ else // density
+ PushDouble( GetGammaDistPDF( fX, fAlpha, fBeta));
}
- else // Verteilung
- PushDouble(GetGammaDist(x, alpha, beta));
}
void ScInterpreter::ScNormInv()
@@ -1575,7 +1728,7 @@ void ScInterpreter::ScGammaInv()
PushInt(0);
else
{
- BOOL bConvError;
+ bool bConvError;
ScGammaDistFunction aFunc( *this, fP, fAlpha, fBeta );
double fStart = fAlpha * fBeta;
double fVal = lcl_IterateInverse( aFunc, fStart*0.5, fStart, bConvError );
@@ -1623,7 +1776,7 @@ void ScInterpreter::ScBetaInv()
PushInt(0);
else
{
- BOOL bConvError;
+ 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 );
@@ -1662,7 +1815,7 @@ void ScInterpreter::ScTInv()
return;
}
- BOOL bConvError;
+ bool bConvError;
ScTDistFunction aFunc( *this, fP, fDF );
double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
if (bConvError)
@@ -1695,7 +1848,7 @@ void ScInterpreter::ScFInv()
return;
}
- BOOL bConvError;
+ bool bConvError;
ScFDistFunction aFunc( *this, fP, fF1, fF2 );
double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
if (bConvError)
@@ -1727,7 +1880,7 @@ void ScInterpreter::ScChiInv()
return;
}
- BOOL bConvError;
+ bool bConvError;
ScChiDistFunction aFunc( *this, fP, fDF );
double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
if (bConvError)