1 /************************************************************************* 2 * 3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. 4 * 5 * Copyright 2000, 2010 Oracle and/or its affiliates. 6 * 7 * OpenOffice.org - a multi-platform office productivity suite 8 * 9 * This file is part of OpenOffice.org. 10 * 11 * OpenOffice.org is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU Lesser General Public License version 3 13 * only, as published by the Free Software Foundation. 14 * 15 * OpenOffice.org is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 * GNU Lesser General Public License version 3 for more details 19 * (a copy is included in the LICENSE file that accompanied this code). 20 * 21 * You should have received a copy of the GNU Lesser General Public License 22 * version 3 along with OpenOffice.org. If not, see 23 * <http://www.openoffice.org/license.html> 24 * for a copy of the LGPLv3 License. 25 * 26 ************************************************************************/ 27 28 // MARKER(update_precomp.py): autogen include statement, do not remove 29 #include "precompiled_sc.hxx" 30 31 // #include <math.h> 32 33 #include <tools/debug.hxx> 34 #include <rtl/logfile.hxx> 35 #include "interpre.hxx" 36 37 double const fHalfMachEps = 0.5 * ::std::numeric_limits<double>::epsilon(); 38 39 // The idea how this group of gamma functions is calculated, is 40 // based on the Cephes library 41 // online http://www.moshier.net/#Cephes [called 2008-02] 42 43 /** You must ensure fA>0.0 && fX>0.0 44 valid results only if fX > fA+1.0 45 uses continued fraction with odd items */ 46 double ScInterpreter::GetGammaContFraction( double fA, double fX ) 47 { 48 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaContFraction" ); 49 50 double const fBigInv = ::std::numeric_limits<double>::epsilon(); 51 double const fBig = 1.0/fBigInv; 52 double fCount = 0.0; 53 double fNum = 0.0; // dummy value 54 double fY = 1.0 - fA; 55 double fDenom = fX + 2.0-fA; 56 double fPk = 0.0; // dummy value 57 double fPkm1 = fX + 1.0; 58 double fPkm2 = 1.0; 59 double fQk = 1.0; // dummy value 60 double fQkm1 = fDenom * fX; 61 double fQkm2 = fX; 62 double fApprox = fPkm1/fQkm1; 63 bool bFinished = false; 64 double fR = 0.0; // dummy value 65 do 66 { 67 fCount = fCount +1.0; 68 fY = fY+ 1.0; 69 fNum = fY * fCount; 70 fDenom = fDenom +2.0; 71 fPk = fPkm1 * fDenom - fPkm2 * fNum; 72 fQk = fQkm1 * fDenom - fQkm2 * fNum; 73 if (fQk != 0.0) 74 { 75 fR = fPk/fQk; 76 bFinished = (fabs( (fApprox - fR)/fR ) <= fHalfMachEps); 77 fApprox = fR; 78 } 79 fPkm2 = fPkm1; 80 fPkm1 = fPk; 81 fQkm2 = fQkm1; 82 fQkm1 = fQk; 83 if (fabs(fPk) > fBig) 84 { 85 // reduce a fraction does not change the value 86 fPkm2 = fPkm2 * fBigInv; 87 fPkm1 = fPkm1 * fBigInv; 88 fQkm2 = fQkm2 * fBigInv; 89 fQkm1 = fQkm1 * fBigInv; 90 } 91 } while (!bFinished && fCount<10000); 92 // most iterations, if fX==fAlpha+1.0; approx sqrt(fAlpha) iterations then 93 if (!bFinished) 94 { 95 SetError(errNoConvergence); 96 } 97 return fApprox; 98 } 99 100 /** You must ensure fA>0.0 && fX>0.0 101 valid results only if fX <= fA+1.0 102 uses power series */ 103 double ScInterpreter::GetGammaSeries( double fA, double fX ) 104 { 105 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaSeries" ); 106 double fDenomfactor = fA; 107 double fSummand = 1.0/fA; 108 double fSum = fSummand; 109 int nCount=1; 110 do 111 { 112 fDenomfactor = fDenomfactor + 1.0; 113 fSummand = fSummand * fX/fDenomfactor; 114 fSum = fSum + fSummand; 115 nCount = nCount+1; 116 } while ( fSummand/fSum > fHalfMachEps && nCount<=10000); 117 // large amount of iterations will be carried out for huge fAlpha, even 118 // if fX <= fAlpha+1.0 119 if (nCount>10000) 120 { 121 SetError(errNoConvergence); 122 } 123 return fSum; 124 } 125 126 /** You must ensure fA>0.0 && fX>0.0) */ 127 double ScInterpreter::GetLowRegIGamma( double fA, double fX ) 128 { 129 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetLowRegIGamma" ); 130 double fLnFactor = fA * log(fX) - fX - GetLogGamma(fA); 131 double fFactor = exp(fLnFactor); // Do we need more accuracy than exp(ln()) has? 132 if (fX>fA+1.0) // includes fX>1.0; 1-GetUpRegIGamma, continued fraction 133 return 1.0 - fFactor * GetGammaContFraction(fA,fX); 134 else // fX<=1.0 || fX<=fA+1.0, series 135 return fFactor * GetGammaSeries(fA,fX); 136 } 137 138 /** You must ensure fA>0.0 && fX>0.0) */ 139 double ScInterpreter::GetUpRegIGamma( double fA, double fX ) 140 { 141 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetUpRegIGamma" ); 142 143 double fLnFactor= fA*log(fX)-fX-GetLogGamma(fA); 144 double fFactor = exp(fLnFactor); //Do I need more accuracy than exp(ln()) has?; 145 if (fX>fA+1.0) // includes fX>1.0 146 return fFactor * GetGammaContFraction(fA,fX); 147 else //fX<=1 || fX<=fA+1, 1-GetLowRegIGamma, series 148 return 1.0 -fFactor * GetGammaSeries(fA,fX); 149 } 150 151 /** Gamma distribution, probability density function. 152 fLambda is "scale" parameter 153 You must ensure fAlpha>0.0 and fLambda>0.0 */ 154 double ScInterpreter::GetGammaDistPDF( double fX, double fAlpha, double fLambda ) 155 { 156 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDistPDF" ); 157 if (fX <= 0.0) 158 return 0.0; // see ODFF 159 else 160 { 161 double fXr = fX / fLambda; 162 // use exp(ln()) only for large arguments because of less accuracy 163 if (fXr > 1.0) 164 { 165 const double fLogDblMax = log( ::std::numeric_limits<double>::max()); 166 if (log(fXr) * (fAlpha-1.0) < fLogDblMax && fAlpha < fMaxGammaArgument) 167 { 168 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha); 169 } 170 else 171 { 172 return exp( (fAlpha-1.0) * log(fXr) - fXr - log(fLambda) - GetLogGamma(fAlpha)); 173 } 174 } 175 else // fXr near to zero 176 { 177 if (fAlpha<fMaxGammaArgument) 178 { 179 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha); 180 } 181 else 182 { 183 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / exp( GetLogGamma(fAlpha)); 184 } 185 } 186 } 187 } 188 189 /** Gamma distribution, cumulative distribution function. 190 fLambda is "scale" parameter 191 You must ensure fAlpha>0.0 and fLambda>0.0 */ 192 double ScInterpreter::GetGammaDist( double fX, double fAlpha, double fLambda ) 193 { 194 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDist" ); 195 if (fX <= 0.0) 196 return 0.0; 197 else 198 return GetLowRegIGamma( fAlpha, fX / fLambda); 199 } 200