Skip to content

C
core
Kaydet (Commit) 5cf55f5b authored tarafından Eike Rathke [er]'s avatar Eike Rathke [er] Kaydeden (comit) Michael Meeks
Dosyalara gözat
Seçenekler

calc66: #i75279# improved accuracy of BINOMDIST; patch from Regina 

Conflicts:
	sc/source/core/tool/interpr3.cxx
üst e21858af
Hide whitespace changes
Inline Side-by-side
Showing with 113 additions and 141 deletions
sc/source/core/inc/interpre.hxx
Dosyayı görüntüle @ 5cf55f5b
...@@ -709,6 +709,7 @@ double GetGamma(double x); ...@@ -709,6 +709,7 @@ double GetGamma(double x);
double GetLogGamma(double x); double GetLogGamma(double x);
double GetBeta(double fAlpha, double fBeta); double GetBeta(double fAlpha, double fBeta);
double GetLogBeta(double fAlpha, double fBeta); double GetLogBeta(double fAlpha, double fBeta);
double GetBinomDistPMF(double x, double n, double p); //probability mass function
void ScLogGamma(); void ScLogGamma();
void ScGamma(); void ScGamma();
void ScPhi(); void ScPhi();
......
sc/source/core/tool/interpr3.cxx
Dosyayı görüntüle @ 5cf55f5b
...@@ -889,17 +889,18 @@ double ScInterpreter::GetBetaDistPDF(double fX, double fA, double fB) ...@@ -889,17 +889,18 @@ double ScInterpreter::GetBetaDistPDF(double fX, double fA, double fB)
const double fLogDblMin = log( ::std::numeric_limits<double>::min()); const double fLogDblMin = log( ::std::numeric_limits<double>::min());
double fLogY = (fX < 0.1) ? ::rtl::math::log1p(-fX) : log(0.5-fX+0.5); double fLogY = (fX < 0.1) ? ::rtl::math::log1p(-fX) : log(0.5-fX+0.5);
double fLogX = log(fX); double fLogX = log(fX);
double fAm1 = fA-1.0; double fAm1LogX = (fA-1.0) * fLogX;
double fBm1 = fB-1.0; double fBm1LogY = (fB-1.0) * fLogY;
double fLogBeta = GetLogBeta(fA,fB); double fLogBeta = GetLogBeta(fA,fB);
// check whether parts over- or underflow // check whether parts over- or underflow
if ( fAm1 * fLogX < fLogDblMax && fAm1 * fLogX > fLogDblMin if ( fAm1LogX < fLogDblMax && fAm1LogX > fLogDblMin
&& fBm1 * fLogY < fLogDblMax && fBm1* fLogY > fLogDblMin && fBm1LogY < fLogDblMax && fBm1LogY > fLogDblMin
&& fLogBeta < fLogDblMax && fLogBeta > 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); return pow(fX,fA-1.0) * pow(0.5-fX+0.5,fB-1.0) / GetBeta(fA,fB);
else // need logarithm; else // need logarithm;
// might overflow as a whole, but seldom, not worth to pre-detect it // might overflow as a whole, but seldom, not worth to pre-detect it
return exp((fA-1.0)*fLogX + (fB-1.0)* fLogY - fLogBeta); return exp( fAm1LogX + fBm1LogY - fLogBeta);
} }
/* /*
...@@ -1165,121 +1166,106 @@ void ScInterpreter::ScVariationen2() ...@@ -1165,121 +1166,106 @@ void ScInterpreter::ScVariationen2()
} }
} }
void ScInterpreter::ScB()
{ double ScInterpreter::GetBinomDistPMF(double x, double n, double p)
RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScB" ); // used in ScB and ScBinomDist
sal_uInt8 nParamCount = GetByte(); // preconditions: 0.0 <= x <= n, 0.0 < p < 1.0; x,n integral although double
if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
return ;
if (nParamCount == 3)
{
double x = ::rtl::math::approxFloor(GetDouble());
double p = GetDouble();
double n = ::rtl::math::approxFloor(GetDouble());
if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
PushIllegalArgument();
else
{ {
double q = 1.0 - p; double q = (0.5 - p) + 0.5;
double fFactor = pow(q, n); double fFactor = pow(q, n);
if (fFactor == 0.0) if (fFactor <=::std::numeric_limits<double>::min())
{ {
fFactor = pow(p, n); fFactor = pow(p, n);
if (fFactor == 0.0) if (fFactor <= ::std::numeric_limits<double>::min())
PushNoValue(); return GetBetaDistPDF(p, x+1.0, n-x+1.0)/(n+1.0);
else else
{ {
sal_uLong max = (sal_uLong) (n - x); sal_uInt32 max = static_cast<sal_uInt32>(n - x);
for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
fFactor *= (n-i)/(i+1)*q/p; fFactor *= (n-i)/(i+1)*q/p;
PushDouble(fFactor); return fFactor;
} }
} }
else else
{ {
sal_uLong max = (sal_uLong) x; sal_uInt32 max = static_cast<sal_uInt32>(x);
for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
fFactor *= (n-i)/(i+1)*p/q; fFactor *= (n-i)/(i+1)*p/q;
PushDouble(fFactor); return fFactor;
}
} }
} }
else if (nParamCount == 4)
{ double lcl_GetBinomDistRange(double n, double xs,double xe,
double xe = GetDouble(); double fFactor /* q^n */, double p, double q)
double xs = GetDouble(); //preconditions: 0.0 <= xs < xe <= n; xs,xe,n integral although double
double p = GetDouble();
double n = GetDouble();
// alter Stand 300-SC
// if ((xs < n) && (xe < n) && (p < 1.0))
// {
// double Varianz = sqrt(n * p * (1.0 - p));
// xs = fabs(xs - (n * p /* / 2.0 STE */ ));
// xe = fabs(xe - (n * p /* / 2.0 STE */ ));
//// STE double nVal = gauss((xs + 0.5) / Varianz) + gauss((xe + 0.5) / Varianz);
// double nVal = fabs(gauss(xs / Varianz) - gauss(xe / Varianz));
// PushDouble(nVal);
// }
bool bIsValidX = ( 0.0 <= xs && xs <= xe && xe <= n);
if ( bIsValidX && 0.0 < p && p < 1.0 )
{
double q = 1.0 - p;
double fFactor = pow(q, n);
if (fFactor == 0.0)
{
fFactor = pow(p, n);
if (fFactor == 0.0)
PushNoValue();
else
{ {
double fSum = 0.0; sal_uInt32 i;
sal_uLong max; double fSum;
if (xe < (sal_uLong) n) // skip summands index 0 to xs-1, start sum with index xs
max = (sal_uLong) (n-xe)-1; sal_uInt32 nXs = static_cast<sal_uInt32>( xs );
else for (i = 1; i <= nXs && fFactor > 0.0; i++)
max = 0; fFactor *= (n-i+1)/i * p/q;
sal_uLong i; fSum = fFactor; // Summand xs
for (i = 0; i < max && fFactor > 0.0; i++) sal_uInt32 nXe = static_cast<sal_uInt32>(xe);
fFactor *= (n-i)/(i+1)*q/p; for (i = nXs+1; i <= nXe && fFactor > 0.0; i++)
if (xs < (sal_uLong) n)
max = (sal_uLong) (n-xs);
else
fSum = fFactor;
for (; i < max && fFactor > 0.0; i++)
{ {
fFactor *= (n-i)/(i+1)*q/p; fFactor *= (n-i+1)/i * p/q;
fSum += fFactor; fSum += fFactor;
} }
PushDouble(fSum); return (fSum>1.0) ? 1.0 : fSum;
} }
void ScInterpreter::ScB()
{
RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScB" );
sal_uInt8 nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
return ;
if (nParamCount == 3) // 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 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)
{ {
sal_uLong max; if (xs == xe) // mass function
double fSum; PushDouble(GetBinomDistPMF(xs,n,p));
if ( (sal_uLong) xs == 0) else
{ {
fSum = fFactor; double fFactor = pow(q, n);
max = 0; if (fFactor > ::std::numeric_limits<double>::min())
} PushDouble(lcl_GetBinomDistRange(n,xs,xe,fFactor,p,q));
else else
{ {
max = (sal_uLong) xs-1; fFactor = pow(p, n);
fSum = 0.0; 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));
} }
sal_uLong i;
for (i = 0; i < max && fFactor > 0.0; i++)
fFactor *= (n-i)/(i+1)*p/q;
if ((sal_uLong)xe == 0) // beide 0
fSum = fFactor;
else else
max = (sal_uLong) xe; PushDouble(GetBetaDist(q,n-xe,xe+1.0)-GetBetaDist(q,n-xs+1,xs) );
for (; i < max && fFactor > 0.0; i++)
{
fFactor *= (n-i)/(i+1)*p/q;
fSum += fFactor;
} }
PushDouble(fSum);
} }
} }
else else
...@@ -1304,78 +1290,63 @@ void ScInterpreter::ScBinomDist() ...@@ -1304,78 +1290,63 @@ void ScInterpreter::ScBinomDist()
RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScBinomDist" ); RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScBinomDist" );
if ( MustHaveParamCount( GetByte(), 4 ) ) if ( MustHaveParamCount( GetByte(), 4 ) )
{ {
double kum = GetDouble(); // 0 oder 1 bool bIsCum = GetBool(); // false=mass function; true=cumulative
double p = GetDouble(); // p double p = GetDouble();
double n = ::rtl::math::approxFloor(GetDouble()); // n double n = ::rtl::math::approxFloor(GetDouble());
double x = ::rtl::math::approxFloor(GetDouble()); // x double x = ::rtl::math::approxFloor(GetDouble());
double fFactor, q; double q = (0.5 - p) + 0.5; // get one bit more for p near 1.0
double fFactor, fSum;
if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0) if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
PushIllegalArgument();
else if (kum == 0.0) // Dichte
{ {
q = 1.0 - p; PushIllegalArgument();
fFactor = pow(q, n); return;
if (fFactor == 0.0)
{
fFactor = pow(p, n);
if (fFactor == 0.0)
PushNoValue();
else
{
sal_uLong max = (sal_uLong) (n - x);
for (sal_uLong i = 0; i < max && fFactor > 0.0; i++)
fFactor *= (n-i)/(i+1)*q/p;
PushDouble(fFactor);
}
} }
else if ( p == 0.0)
{ {
sal_uLong max = (sal_uLong) x; PushDouble( (x==0.0 || bIsCum) ? 1.0 : 0.0 );
for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) return;
fFactor *= (n-i)/(i+1)*p/q;
PushDouble(fFactor);
} }
if ( p == 1.0)
{
PushDouble( (x==n) ? 1.0 : 0.0);
return;
} }
else // Verteilung if (!bIsCum)
PushDouble( GetBinomDistPMF(x,n,p));
else
{ {
if (n == x) if (x == n)
PushDouble(1.0); PushDouble(1.0);
else else
{ {
double fSum;
q = 1.0 - p;
fFactor = pow(q, n); fFactor = pow(q, n);
if (fFactor == 0.0) if (x == 0.0)
PushDouble(fFactor);
else
if (fFactor <= ::std::numeric_limits<double>::min())
{ {
fFactor = pow(p, n); fFactor = pow(p, n);
if (fFactor == 0.0) if (fFactor <= ::std::numeric_limits<double>::min())
PushNoValue(); PushDouble(GetBetaDist(q,n-x,x+1.0));
else else
{ {
if (fFactor > fMachEps)
{
fSum = 1.0 - fFactor; fSum = 1.0 - fFactor;
sal_uLong max = (sal_uLong) (n - x) - 1; sal_uInt32 max = static_cast<sal_uInt32> (n - x) - 1;
for (sal_uLong i = 0; i < max && fFactor > 0.0; i++) for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
{ {
fFactor *= (n-i)/(i+1)*q/p; fFactor *= (n-i)/(i+1)*q/p;
fSum -= fFactor; fSum -= fFactor;
} }
if (fSum < 0.0) PushDouble( (fSum < 0.0) ? 0.0 : fSum );
PushDouble(0.0);
else
PushDouble(fSum);
}
} }
else else
{ PushDouble(lcl_GetBinomDistRange(n,n-x,n,fFactor,q,p));
fSum = fFactor;
sal_uLong max = (sal_uLong) x;
for (sal_uLong i = 0; i < max && fFactor > 0.0; i++)
{
fFactor *= (n-i)/(i+1)*p/q;
fSum += fFactor;
} }
PushDouble(fSum);
} }
else
PushDouble( lcl_GetBinomDistRange(n,0.0,x,fFactor,p,q)) ;
} }
} }
} }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment