Branch data Line data Source code
1 : : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : : /*************************************************************************
3 : : *
4 : : * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
5 : : *
6 : : * Copyright 2000, 2010 Oracle and/or its affiliates.
7 : : *
8 : : * OpenOffice.org - a multi-platform office productivity suite
9 : : *
10 : : * This file is part of OpenOffice.org.
11 : : *
12 : : * OpenOffice.org is free software: you can redistribute it and/or modify
13 : : * it under the terms of the GNU Lesser General Public License version 3
14 : : * only, as published by the Free Software Foundation.
15 : : *
16 : : * OpenOffice.org is distributed in the hope that it will be useful,
17 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 : : * GNU Lesser General Public License version 3 for more details
20 : : * (a copy is included in the LICENSE file that accompanied this code).
21 : : *
22 : : * You should have received a copy of the GNU Lesser General Public License
23 : : * version 3 along with OpenOffice.org. If not, see
24 : : * <http://www.openoffice.org/license.html>
25 : : * for a copy of the LGPLv3 License.
26 : : *
27 : : ************************************************************************/
28 : :
29 : :
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 : 0 : double ScInterpreter::GetGammaContFraction( double fA, double fX )
43 : : {
44 : : RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaContFraction" );
45 : :
46 : 0 : double const fBigInv = ::std::numeric_limits<double>::epsilon();
47 : 0 : double const fBig = 1.0/fBigInv;
48 : 0 : double fCount = 0.0;
49 : 0 : double fNum = 0.0; // dummy value
50 : 0 : double fY = 1.0 - fA;
51 : 0 : double fDenom = fX + 2.0-fA;
52 : 0 : double fPk = 0.0; // dummy value
53 : 0 : double fPkm1 = fX + 1.0;
54 : 0 : double fPkm2 = 1.0;
55 : 0 : double fQk = 1.0; // dummy value
56 : 0 : double fQkm1 = fDenom * fX;
57 : 0 : double fQkm2 = fX;
58 : 0 : double fApprox = fPkm1/fQkm1;
59 : 0 : bool bFinished = false;
60 : 0 : double fR = 0.0; // dummy value
61 [ # # ][ # # ]: 0 : do
[ # # ]
62 : : {
63 : 0 : fCount = fCount +1.0;
64 : 0 : fY = fY+ 1.0;
65 : 0 : fNum = fY * fCount;
66 : 0 : fDenom = fDenom +2.0;
67 : 0 : fPk = fPkm1 * fDenom - fPkm2 * fNum;
68 : 0 : fQk = fQkm1 * fDenom - fQkm2 * fNum;
69 [ # # ]: 0 : if (fQk != 0.0)
70 : : {
71 : 0 : fR = fPk/fQk;
72 : 0 : bFinished = (fabs( (fApprox - fR)/fR ) <= fHalfMachEps);
73 : 0 : fApprox = fR;
74 : : }
75 : 0 : fPkm2 = fPkm1;
76 : 0 : fPkm1 = fPk;
77 : 0 : fQkm2 = fQkm1;
78 : 0 : fQkm1 = fQk;
79 [ # # ]: 0 : if (fabs(fPk) > fBig)
80 : : {
81 : : // reduce a fraction does not change the value
82 : 0 : fPkm2 = fPkm2 * fBigInv;
83 : 0 : fPkm1 = fPkm1 * fBigInv;
84 : 0 : fQkm2 = fQkm2 * fBigInv;
85 : 0 : fQkm1 = fQkm1 * fBigInv;
86 : : }
87 : 0 : } while (!bFinished && fCount<10000);
88 : : // most iterations, if fX==fAlpha+1.0; approx sqrt(fAlpha) iterations then
89 [ # # ]: 0 : if (!bFinished)
90 : : {
91 : 0 : SetError(errNoConvergence);
92 : : }
93 : 0 : 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 : 0 : double ScInterpreter::GetGammaSeries( double fA, double fX )
100 : : {
101 : : RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaSeries" );
102 : 0 : double fDenomfactor = fA;
103 : 0 : double fSummand = 1.0/fA;
104 : 0 : double fSum = fSummand;
105 : 0 : int nCount=1;
106 [ # # ][ # # ]: 0 : do
[ # # ]
107 : : {
108 : 0 : fDenomfactor = fDenomfactor + 1.0;
109 : 0 : fSummand = fSummand * fX/fDenomfactor;
110 : 0 : fSum = fSum + fSummand;
111 : 0 : 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 [ # # ]: 0 : if (nCount>10000)
116 : : {
117 : 0 : SetError(errNoConvergence);
118 : : }
119 : 0 : return fSum;
120 : : }
121 : :
122 : : /** You must ensure fA>0.0 && fX>0.0) */
123 : 0 : double ScInterpreter::GetLowRegIGamma( double fA, double fX )
124 : : {
125 : : RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetLowRegIGamma" );
126 : 0 : double fLnFactor = fA * log(fX) - fX - GetLogGamma(fA);
127 : 0 : double fFactor = exp(fLnFactor); // Do we need more accuracy than exp(ln()) has?
128 [ # # ]: 0 : if (fX>fA+1.0) // includes fX>1.0; 1-GetUpRegIGamma, continued fraction
129 : 0 : return 1.0 - fFactor * GetGammaContFraction(fA,fX);
130 : : else // fX<=1.0 || fX<=fA+1.0, series
131 : 0 : return fFactor * GetGammaSeries(fA,fX);
132 : : }
133 : :
134 : : /** You must ensure fA>0.0 && fX>0.0) */
135 : 0 : double ScInterpreter::GetUpRegIGamma( double fA, double fX )
136 : : {
137 : : RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetUpRegIGamma" );
138 : :
139 : 0 : double fLnFactor= fA*log(fX)-fX-GetLogGamma(fA);
140 : 0 : double fFactor = exp(fLnFactor); //Do I need more accuracy than exp(ln()) has?;
141 [ # # ]: 0 : if (fX>fA+1.0) // includes fX>1.0
142 : 0 : return fFactor * GetGammaContFraction(fA,fX);
143 : : else //fX<=1 || fX<=fA+1, 1-GetLowRegIGamma, series
144 : 0 : 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 : 0 : double ScInterpreter::GetGammaDistPDF( double fX, double fAlpha, double fLambda )
151 : : {
152 : : RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDistPDF" );
153 [ # # ]: 0 : if (fX < 0.0)
154 : 0 : return 0.0; // see ODFF
155 [ # # ]: 0 : else if (fX == 0)
156 : : // in this case 0^0 isn't zero
157 : : {
158 [ # # ]: 0 : if (fAlpha < 1.0)
159 : : {
160 : 0 : SetError(errDivisionByZero); // should be #DIV/0
161 : 0 : return HUGE_VAL;
162 : : }
163 [ # # ]: 0 : else if (fAlpha == 1)
164 : : {
165 : 0 : return (1.0 / fLambda);
166 : : }
167 : : else
168 : : {
169 : 0 : return 0.0;
170 : : }
171 : : }
172 : : else
173 : : {
174 : 0 : double fXr = fX / fLambda;
175 : : // use exp(ln()) only for large arguments because of less accuracy
176 [ # # ]: 0 : if (fXr > 1.0)
177 : : {
178 : 0 : const double fLogDblMax = log( ::std::numeric_limits<double>::max());
179 [ # # ][ # # ]: 0 : if (log(fXr) * (fAlpha-1.0) < fLogDblMax && fAlpha < fMaxGammaArgument)
[ # # ]
180 : : {
181 : 0 : return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha);
182 : : }
183 : : else
184 : : {
185 : 0 : return exp( (fAlpha-1.0) * log(fXr) - fXr - log(fLambda) - GetLogGamma(fAlpha));
186 : : }
187 : : }
188 : : else // fXr near to zero
189 : : {
190 [ # # ]: 0 : if (fAlpha<fMaxGammaArgument)
191 : : {
192 : 0 : return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / GetGamma(fAlpha);
193 : : }
194 : : else
195 : : {
196 : 0 : return pow( fXr, fAlpha-1.0) * exp(-fXr) / fLambda / exp( GetLogGamma(fAlpha));
197 : : }
198 : : }
199 : : }
200 : : }
201 : :
202 : : /** Gamma distribution, cumulative distribution function.
203 : : fLambda is "scale" parameter
204 : : You must ensure fAlpha>0.0 and fLambda>0.0 */
205 : 0 : double ScInterpreter::GetGammaDist( double fX, double fAlpha, double fLambda )
206 : : {
207 : : RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetGammaDist" );
208 [ # # ]: 0 : if (fX <= 0.0)
209 : 0 : return 0.0;
210 : : else
211 : 0 : return GetLowRegIGamma( fAlpha, fX / fLambda);
212 : : }
213 : :
214 : : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
