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