1 /************************************************************** 2 * 3 * Licensed to the Apache Software Foundation (ASF) under one 4 * or more contributor license agreements. See the NOTICE file 5 * distributed with this work for additional information 6 * regarding copyright ownership. The ASF licenses this file 7 * to you under the Apache License, Version 2.0 (the 8 * "License"); you may not use this file except in compliance 9 * with the License. You may obtain a copy of the License at 10 * 11 * http://www.apache.org/licenses/LICENSE-2.0 12 * 13 * Unless required by applicable law or agreed to in writing, 14 * software distributed under the License is distributed on an 15 * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY 16 * KIND, either express or implied. See the License for the 17 * specific language governing permissions and limitations 18 * under the License. 19 * 20 *************************************************************/ 21 22 23 24 // MARKER(update_precomp.py): autogen include statement, do not remove 25 #include "precompiled_sc.hxx" 26 27 // #include <math.h> 28 29 #include <tools/debug.hxx> 30 #include <rtl/logfile.hxx> 31 #include "interpre.hxx" 32 33 double const fHalfMachEps = 0.5 * ::std::numeric_limits<double>::epsilon(); 34 35 // The idea how this group of gamma functions is calculated, is 36 // based on the Cephes library 37 // online http://www.moshier.net/#Cephes [called 2008-02] 38 39 /** You must ensure fA>0.0 && fX>0.0 40 valid results only if fX > fA+1.0 41 uses continued fraction with odd items */ 42 double ScInterpreter::GetGammaContFraction( double fA, double fX ) 43 { 44 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaContFraction" ); 45 46 double const fBigInv = ::std::numeric_limits<double>::epsilon(); 47 double const fBig = 1.0/fBigInv; 48 double fCount = 0.0; 49 double fNum = 0.0; // dummy value 50 double fY = 1.0 - fA; 51 double fDenom = fX + 2.0-fA; 52 double fPk = 0.0; // dummy value 53 double fPkm1 = fX + 1.0; 54 double fPkm2 = 1.0; 55 double fQk = 1.0; // dummy value 56 double fQkm1 = fDenom * fX; 57 double fQkm2 = fX; 58 double fApprox = fPkm1/fQkm1; 59 bool bFinished = false; 60 double fR = 0.0; // dummy value 61 do 62 { 63 fCount = fCount +1.0; 64 fY = fY+ 1.0; 65 fNum = fY * fCount; 66 fDenom = fDenom +2.0; 67 fPk = fPkm1 * fDenom - fPkm2 * fNum; 68 fQk = fQkm1 * fDenom - fQkm2 * fNum; 69 if (fQk != 0.0) 70 { 71 fR = fPk/fQk; 72 bFinished = (fabs( (fApprox - fR)/fR ) <= fHalfMachEps); 73 fApprox = fR; 74 } 75 fPkm2 = fPkm1; 76 fPkm1 = fPk; 77 fQkm2 = fQkm1; 78 fQkm1 = fQk; 79 if (fabs(fPk) > fBig) 80 { 81 // reduce a fraction does not change the value 82 fPkm2 = fPkm2 * fBigInv; 83 fPkm1 = fPkm1 * fBigInv; 84 fQkm2 = fQkm2 * fBigInv; 85 fQkm1 = fQkm1 * fBigInv; 86 } 87 } while (!bFinished && fCount<10000); 88 // most iterations, if fX==fAlpha+1.0; approx sqrt(fAlpha) iterations then 89 if (!bFinished) 90 { 91 SetError(errNoConvergence); 92 } 93 return fApprox; 94 } 95 96 /** You must ensure fA>0.0 && fX>0.0 97 valid results only if fX <= fA+1.0 98 uses power series */ 99 double ScInterpreter::GetGammaSeries( double fA, double fX ) 100 { 101 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaSeries" ); 102 double fDenomfactor = fA; 103 double fSummand = 1.0/fA; 104 double fSum = fSummand; 105 int nCount=1; 106 do 107 { 108 fDenomfactor = fDenomfactor + 1.0; 109 fSummand = fSummand * fX/fDenomfactor; 110 fSum = fSum + fSummand; 111 nCount = nCount+1; 112 } while ( fSummand/fSum > fHalfMachEps && nCount<=10000); 113 // large amount of iterations will be carried out for huge fAlpha, even 114 // if fX <= fAlpha+1.0 115 if (nCount>10000) 116 { 117 SetError(errNoConvergence); 118 } 119 return fSum; 120 } 121 122 /** You must ensure fA>0.0 && fX>0.0) */ 123 double ScInterpreter::GetLowRegIGamma( double fA, double fX ) 124 { 125 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetLowRegIGamma" ); 126 double fLnFactor = fA * log(fX) - fX - GetLogGamma(fA); 127 double fFactor = exp(fLnFactor); // Do we need more accuracy than exp(ln()) has? 128 if (fX>fA+1.0) // includes fX>1.0; 1-GetUpRegIGamma, continued fraction 129 return 1.0 - fFactor * GetGammaContFraction(fA,fX); 130 else // fX<=1.0 || fX<=fA+1.0, series 131 return fFactor * GetGammaSeries(fA,fX); 132 } 133 134 /** You must ensure fA>0.0 && fX>0.0) */ 135 double ScInterpreter::GetUpRegIGamma( double fA, double fX ) 136 { 137 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetUpRegIGamma" ); 138 139 double fLnFactor= fA*log(fX)-fX-GetLogGamma(fA); 140 double fFactor = exp(fLnFactor); //Do I need more accuracy than exp(ln()) has?; 141 if (fX>fA+1.0) // includes fX>1.0 142 return fFactor * GetGammaContFraction(fA,fX); 143 else //fX<=1 || fX<=fA+1, 1-GetLowRegIGamma, series 144 return 1.0 -fFactor * GetGammaSeries(fA,fX); 145 } 146 147 /** Gamma distribution, probability density function. 148 fLambda is "scale" parameter 149 You must ensure fAlpha>0.0 and fLambda>0.0 */ 150 double ScInterpreter::GetGammaDistPDF( double fX, double fAlpha, double fLambda ) 151 { 152 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDistPDF" ); 153 if (fX <= 0.0) 154 return 0.0; // see ODFF 155 else 156 { 157 double fXr = fX / fLambda; 158 // use exp(ln()) only for large arguments because of less accuracy 159 if (fXr > 1.0) 160 { 161 const double fLogDblMax = log( ::std::numeric_limits<double>::max()); 162 if (log(fXr) * (fAlpha-1.0) < fLogDblMax && fAlpha < fMaxGammaArgument) 163 { 164 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha); 165 } 166 else 167 { 168 return exp( (fAlpha-1.0) * log(fXr) - fXr - log(fLambda) - GetLogGamma(fAlpha)); 169 } 170 } 171 else // fXr near to zero 172 { 173 if (fAlpha<fMaxGammaArgument) 174 { 175 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha); 176 } 177 else 178 { 179 return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / exp( GetLogGamma(fAlpha)); 180 } 181 } 182 } 183 } 184 185 /** Gamma distribution, cumulative distribution function. 186 fLambda is "scale" parameter 187 You must ensure fAlpha>0.0 and fLambda>0.0 */ 188 double ScInterpreter::GetGammaDist( double fX, double fAlpha, double fLambda ) 189 { 190 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDist" ); 191 if (fX <= 0.0) 192 return 0.0; 193 else 194 return GetLowRegIGamma( fAlpha, fX / fLambda); 195 } 196