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 */
GetGammaContFraction(double fA,double fX)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 */
GetGammaSeries(double fA,double fX)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) */
GetLowRegIGamma(double fA,double fX)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) */
GetUpRegIGamma(double fA,double fX)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 */
GetGammaDistPDF(double fX,double fAlpha,double fLambda)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 */
GetGammaDist(double fX,double fAlpha,double fLambda)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