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